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Abstract 

Macroscale  rate-limited  sorption  modeling  was  tested  using  a  production  transport  code, 
the  GMS/FEMW\TER  ground-water  modeling  package.  The  code  (\fersion  1.1  ofFEMWATER. 
dated  1  August  1 995)  was  applied  to  a  3D  conceptual  model  developed  from  a  field  site  at  Dover 
AFB,  DL.  A  simulation  was  performed  of  a  200  hour  contaminant  injection  pulse  followed  by 
clean  water  flushing.  A  moment  analysis  performed  on  the  resulting  breakthrough  curve  vali¬ 
dated  code  self-consistency.  Another  injection  pulse  simulation  showed  that  retardation  tempo¬ 
rally  delays  the  breakthrough  peak.  Transport  simulations  of  pulsed  clean  water  pumping  of  the 
test  cell  with  a  prescribed  initial  contaminant  distribution  demonstrated  both  tailing  and  rebound 
without  any  additional  microscale  modeling.  In  comparison  with  both  previous  numerical  so¬ 
lutions  and  the  actual  field  data  from  the  Dover  AFB  test  site,  FEMWATER  has  demonstrated 
high  numerical  dispersivity.  For  an  initial  contaminant  distribution  corresponding  to  the  field 
data,  the  FEMWATER  breakthrough  curve  was  much  flatter  than  the  experimental  result,  failing 
to  capture  the  plug-like  elution  of  the  field  site. 
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THREE-DIMENSIONAL  CALCULATION  OF 
CONTAMINANT  TRANSPORT  IN  GROUNDWATER  AT 

A  DOVER  AFB  SITE 


1.  INTRODUCTION 


1.1  BACKGROUND 

l.l.a  THE  ENVIRONMENT 

We  live  in  an  era  of  ever-increasing  efforts  to  right  the  environmental  abuses  of  yesteryear. 
A  major  area  of  concern  is  the  condition  of  our  ground  water -the  entity  that  suffers  when  offen¬ 
sive  substances  are  dumped  ”below  decks”[30,  pg  463]  .  According  to  the  1989  Toxic  Release 
Inventory  by  the  Environmental  Protection  Agency  (ER\),  14%  of  toxic  chemical  releases  were 
underground  [10,  pg  173] .  Almost  all  rural  households  (95%)  and  50%  of  the  general  popula¬ 
tion  depend  upon  ground  water  as  their  primary  drinking  water  source  [23,  pg  250]  .  Ground- 
water  contamination  exists  at  more  than  85%  of  the  1 208  sites  included  on  the  National  Priority 
List  (NPL).  Further,  over  33,000  other  sites  have  been  identified  and  included  in  the  Compre¬ 
hensive  Environmental  Response,  Compensation,  and  Liability  Information  System  for  ranking 
and  potential  inclusion  on  the  NPL.  The  EPA  has  identified  or  suspected  contaminant  release  to 
ground  water  at  more  than  1700  Resource  Conservation  and  Recovery  Act  facilities.  Against 
this  backdrop  comes  the  disturbing  realization  that  the  Air  Force  has  been  a  prime  polluting  cul¬ 
prit  in  many  cases.  Typical  scenarios  include  leaky  underground  fuel  storage  tanks  and  careless 
disposal  of  toxic  industrial  chemicals  like  solvents  and  paint  strippers.  The  Air  Force  is  engaged 
in  a  program  to  identify,  assess,  and  remediate  hazardous  waste  sites  at  Air  Force  installations 
throughout  the  United  States[12,  pg  24]  .  These  efforts  are  driven  partially  by  the  reality  of 
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impending  base  closures  and  the  subsequent  turn-over  of  potentially  contaminated  lands  to  the 
civilian  sector. 

l.l.b  PUMP-AND-TREAT 

”Pump-and-Treat”  (PAT)  is  the  collective  term  used  to  describe  contaminated  site  reme¬ 
diation  schemes  in  which  water  is  extracted  from  the  aquifer  by  wells  or  drains  followed  by 
treatment  of  the  extracted  water  [25,  pg  630] ,  [12,  pg  25] ,  [18,  pg  119] ,  [19,  pg  216]  .  PAT 
remediation  is  the  most  frequently  used  remediation  method  used  in  practice  [33,  pg  1464] . 

Field  experience  with  PAT  has  demonstrated  several  trends  in  the  effectiveness  of  the 
process  [25,  pg  630] ,  [18,  pg  119] ,  [29,  pg  44] : 

1 .  Containment  of  ground-water  plumes  was  usually  achieved; 

2.  Contaminant  concentrations  dropped  significantly,  initially,  followed  by  a  leveling  out  and; 

3.  After  the  period  of  initial  decline,  the  continued  decreases  in  concentration  were  usually 

slower  than  anticipated  (sometimes  by  decades).  This  effect  is  known  as  "tailing”. 

4.  After  the  cessation  of  pumping,  there  were  local  increases  in  contaminant  concentration. 

An  effect  known  as  ’’rebound”. 

Tailing  and  rebound  are  both  attributed  to  "rate-limited  sorption”,  also  called  ’’nonequi¬ 
librium  sorption”.  The  term  sorption  describes  processes  that  exchange  solutes  between  the 
fluid  and  solid  phases  of  a  medium.  It  is  the  generic  term  used  to  encompass  the  phenomena 
of  adsorption,  desorption,  and  other  processes.  Equilibrium  assumes  that  the  time  scale  of  the 
chemical  binding  and  release  of  solute  onto  porous  and  solid  particles  within  the  soil  matrix 
is  much  faster  than  the  time  scale  of  bulk  flow  advection.  So  equilibrium  sorption  means  that 
solute  moves  to  or  from  the  solid  matrix  as  quickly  as  it  advects  downstream.  Non-equilibrium 
sorption  arises  when  physical  or  chemical  processes  at  the  single  pore  level  are  slow  relative 
to  advection  in  the  bulk  media  [32,  pg  1]  .  The  rate  of  sorption  is  controlled  by  one  of  three 
possible  rate-limiting  steps  [32,  pg  9] : 
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1.  chemical  binding  and  release, 

2.  diffusion  through  immobile  fluid,  and, 

3.  diffusion  through  organic  matter 

According  to  several  recent  studies  [35,  pg  1699] ,  non-equilibrium  solute  transport  is  a 
key  limitation  to  contaminant  removal  from  ground  water  by  ’’pump  and  treat”  methods.  Slow 
sorption  (with  attendant  tailing  and  rebound  effects)  of  organic  substances  within  the  soil  has 
resulted  in  long  periods  of  treatment,  almost  always  in  excess  of  those  predicted  by  conventional 
equilibrium  modeling,  with  attendant  costs.  Additionally,  such  scenarios  lead  to  potentially 
difficult  pumping  strategy  determinations.  As  a  result,  continuous  pumping  did  not  prove  to  be 
practical  in  many  situations  [25,  pg  630] ,  [18,  pg  119] ,  [29,  pg  44]  .  ’’Pulse  pumping”  refers 
to  intermittent  pumping  processes  designed  to  improve  the  plume  capture  and  reduce  quantities 
of  water  extracted,  to  reach  the  goal  of ’’clean  water”  more  efficiently.  Field  experiments  and 
modeling  studies  have  demonstrated  that  this  optimization  of  the  process  may  achieve  favorable 
results  [3,  pg  165] ,  [18,  pg  122] ,  [12,  pg25] ,  [17,  pg  37] ,  [1,  pg  53] .  Hence  the  importance  of 
accurate  numerical  modeling  of  the  ground-water  contaminant  transport  process:  remediation 
engineers  can  rely  on  these  models  to  help  optimize  their  pumping  schedules  -  saving  money 
in  the  long  run. 

l.l.c  MODELING  CONCERNS 

As  in  any  computational  effort,  the  emphasis  is  on  gleaning  as  much  physical  insight  from 
the  model  of  a  process  or  system,  while  using  the  least  computer  resources  possible,  thus  min¬ 
imizing  cost.  Often,  one  cannot  obtain  precise  quantitative  data  but  can  nevertheless  surmise 
important  trends  in  the  way  the  system  responds  to  variations  in  any  inherent  parameters.  More¬ 
over,  the  ability  to  run  numerical  simulations  then  gives  the  researcher  the  platform  to  test  his 
or  her  notions  about  the  system  and  thus  systematically  arrive  at  a  far  narrower  set  of  optimum 
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solutions  than  would  otherwise  be  possible  through  mere  contemplation  or  observation  of  the 
physical  system. 

Having  stated  the  benefits  as  well  as  limits  of  numerical  models,  we  can  now  state  a  key 
requirement  concerning  the  credibility  of  the  model:  the  need  for  validation.  Validation  is  the 
process  by  which  a  researcher  establishes  a  numerical  code’s  credibility  through  a  systematic 
procedure  of  comparing  results  produced  by  the  code  to  those  from  other  sources.  There  are  at 
least  two  types  of  validation  the  researcher  should  perform: 

1 .  Some  sort  of  cross-check  to  determine  whether  the  model  is  behaving  in  a  self-consistent 
manner 

2.  Some  type  of  comparison  to  actual  physical  data.  Here,  governing  parameters  should  match 
one-to-one  between  the  physical  case  being  simulated  and  the  numerical  model’s  inputs. 

1 .1  .d  BRIEF  OVERVIEW  OF  PREVIOUS  ANALYSES 

In  coarse-grained,  homogeneous  aquifers,  advective  transport  dominates  the  transport  process. 
For  this  reason  as  well  as  for  the  simplification  of  the  mathematical  model  in  the  design  of  cap¬ 
ture  and  containment  systems,  it  has  been  common  to  treat  advection  as  the  sole  mechanism  for 
contaminant  transport  [15,  pg  42]  .  The  vast  majority  of  ground-water  contaminant  transport 
models  are  based  on  the  assumption  of  instantaneous  sorption  and  desorption  between  the  liq¬ 
uid  and  the  solid  phases.  This  assumption  is  commonly  called  the  local  equilibrium  assumption 
(LEA).  It’s  validity  and  applicability  have  been  documented  in  the  literature  [34] ,  [40] ,  [13] , 
[35] ,  [24] ,  [32,  pg  9] ,  [6,  pg  33-99] ,  [4] ,  [12,  pg  24-25] ,  [39,  pg  499-528] ,  [5,  pg  353-368] 

.  It  is  important  to  realize  that  when  discussing  the  LEA  versus  rate-limited  sorption,  one  must 
identify  two  distinct  sources  of  non-equilibrium,  namely:  1 )  Physical  non-equilibrium,  in  which 
the  overall  sorption  rate  is  controlled  by  the  rate  at  which  the  solute  is  transported  to  and  from 
the  reacting  soil  surfaces  and  2)  Chemical  non-equilibrium,  in  which  the  overall  sorption  rate  is 
equal  to  the  rate  of  reaction  at  the  soil-solution  interfaces.  Clearly,  the  slowest  process  sets  the 
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rate-limiting  standard.  In  the  present  study,  we  consider  only  physical  non-equilibrium.  Exam¬ 
ples  of  chemical  non-equilibrium  modeling  can  be  found  in  [35,  pg  1696] ,  [16,  pg  602-604] . 
In  the  realm  of  physical  models,  considerable  work  has  been  done  recently  to  develop  models 
which  accurately  describe  rate-limited  sorption  effects  [22] ,  [4] ,  [1] ,  [7] ,  [31,  pg  1457-1470] . 
Most  of  these  models  utilize  a  two-zone  description  of  a  porous  medium.  They  implement  non¬ 
equilibrium  sorption  in  various  ways,  but  all  have  considered  non-equilibrium  effects  occurring 
throughout  the  porous  medium  on  the  microscopic  scale.  Huso  [22]  ,  using  finite  elements, 
and  Adams  and  Virmontes  [1,  pg  1-2]  using  an  analytical  approach,  modeled  physical,  non¬ 
equilibrium  in  radial,  pulsed  pumping  tests.  A  notable  recent  work  is  that  of  Caspers  [7]  who 
modified  the  widely  used  USGS  SUTRA  code  [37]  to  incorporate  rate-limited  as  well  as  equi¬ 
librium  sorption  effects.  In  his  work,  he  describes  rate-limitation  by  either  a  first-order  law,  or 
by  Fickian  diffusion  of  contaminant  through  a  spherical  immobile  region.  Casper’s  equilibrium 
methods  under-predicted  rebound,  while  his  first-order  diffusion  simulations  both  under  and 
over-predicted  rebound  within  the  matrix  for  certain  regions  and  were  somewhat  equivalent  to 
Fickian  diffusion  in  equilibrium  regimes  for  cleanup  time  prediction.  Herman’s  work  [21]  in¬ 
vestigated  non-equilibrium  sorption  effects  from  only  specific  zones,  not  throughout  the  porous 
medium  on  the  microscopic  scale  like  Casper’s  simulations.  In  Herman’s  study,  the  USGS  SU¬ 
TRA  code  was  coupled  with  a  two  dimensional  diffusion  code  to  model  non-equilibrium  effects 
from  specific  layers  only.  This  allowed  the  study  of  non-equilibrium  sorption  from  a  macro¬ 
scopic  view  of  specific  zones.  Herman  simulated  PAT  remediation  at  a  field  site  at  Dover  AFB, 
Delaware  and  studied  the  effects  of  pulsed  and  continuous  pumping  within  the  time  frame  of 
the  actual  field  experiment.  Since  both  Casper’s  and  Herman’s  works  were  based  on  SUTRA, 
they  were  both  hamstrung  by  a  fundamental  limitation  of  that  code,  namely,  the  inability  to  vary 
retardation  coefficient  per  heterogeneities  in  the  soil  strata  under  consideration.  As  Herman  re- 
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ports,  the  SUTRA  code  did  not  allow  specification  of  retardation  coefficients  in  his  loam  layers. 
Accordingly,  in  his  model,  the  distribution  constant  was  one-third  what  it  should  have  been  in 
orange  silty  clay  loam  and  two  orders  of  magnitude  less  than  what  it  should  have  been  for  black 
silty  loam.  This  inflexibility  resulted  in  an  unrealistically  high  rate  of  mass  mobilization  from 
the  clay  layers  of  his  contaminated  strata.  Even  though  the  hydraulic  conductivity  was  low  in 
the  clay  layers,  mass  was  still  mobilized  out  of  the  loam  layers  because  the  significant  retarda¬ 
tion  in  the  loam  layers  was  not  accounted  for  in  the  SUTRA  model. 

Both  Caspers  and  Herman  used  the  so-called  "Split  Operator  Approach”  to  make  the  nu¬ 
merical  problem  more  tractable.  This  approach,  detailed  by  Miller  and  Rabideau  [28]  separates 
the  governing  equations  into  transport  (for  the  mobile  zone)  and  reaction  operators  (diifusion 
equations  for  the  immobile  zone)  and  solves  them  sequentially.  This  splitting  allows  the  separa¬ 
tion  of  the  short  time-scale  process  (mobile  zone)  from  the  long  time-scale  process  (immobile). 
Operator  splitting  leads  to  smaller  systems  of  equations,  which  can  be  solved  faster  than  the 
original  coupled  equation. 

l.l.e  MOTimTION  FOR  CURRENT  RESEARCH 

The  three  motivations  for  conducting  the  present  research  will  be  detailed  in  the  following 
sections.  They  are: 

1 .  Macroscale  approach  to  rate-limited  sorption 

2.  Three-dimensional  vs  two-dimensional  modeling 

3.  The  simulation  of  an  actual  field  experiment  being  conducted  at  Dover  AFB  with  particular 
attention  paid  to: 

A.  Pump  and  Treat 

B.  Tailing 

C.  Rebound 
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l.l.e.  1  MACROSCALE  APPROACH  TO  RATE-LIMITED  SORPTION 


Numerical  modelers  want  to  know  if  they  can  conduct  preliminary  investigations  of  actual 
field  conditions  without  having  to  resort  to  complex  microscale  models  which  require  detailed 
knowledge  of  characteristics  like  pore  geometries,  particle  size  distribution,  or  the  split-operator 
approach  in  which  a  diffusive  sub-model  must  be  specially  prepared.  If  the  bulk  of  retardation 
occurs  in  macroscale  soil  structures  (like  layers  or  lenses),  then  a  numerical  formulation  that 
assumes  local  chemical  equilibrium  and  models  diffusion  as  well  as  advection,  while  providing 
the  capability  to  vary  retardation  coefficient  throughout  the  modeled  configuration  should  be 
able  to  simulate  processes  like  rate-limited  sorption  that  are  diffusion  dominated.  This  is  the  so- 
called  macroscale  approach.  The  fundamental  question  then  is  whether,  devoid  of  any  additional 
diffusive  model,  the  current  methodology  is  capable  of  reproducing  the  diffusion-dominated 
phenomena  of  tailing  and  rebound? 

LLe.2  THREE-DIMENSIONAL  MODELING 

It  is  always  true  that  efficiency  dictates  the  need  for  the  cheapest  and  fastest  simulation 
possible,  hence  the  need  for  2D  models  of  3D  configurations.  However,  it  is  also  true  that  some 
geometries  are  strongly  three-dimensional  and  can  only  be  properly  represented  by  a  ”full-up” 
3D  model.  In  light  of  this,  the  analyst  must  have  a  clear  knowledge  of  when  a  2D  model  is 
sufficient,  when  a  3D  model  is  necessary,  and  what  types  of  characteristic  parameters  or  physical 
phenomena  are  the  determining  factors.  Unfortunately,  the  only  way  to  build  this  knowledge  is 
by  developing  a  mosaic  of  2D  and  3D  simulations  of  sufficient  breadth  that  unmistakable  trends 
can  be  discerned,  recorded  and  exploited  by  researchers.  Obviously,  the  greater  the  number  of 
simulations,  the  greater  the  ’’corporate  knowledge”  of  when  one  can  expect  to  get  useful  results 
from  a  2D  model,  and  when  one  must  resort  to  a  full  3D  model. 
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l.l.e.3  DOVER  AFB  FIELD  EXPERIMENT 


One  of  the  objectives  of  this  work  was  to  support  a  pump-and-treat  field  experiment  being 
conducted  at  Dover  AFB.  This  site  was  specifically  prepared  to  test  pulse-pumping  etficiency 
with  particular  attention  paid  to  tailing  and  rebound.  It’s  soil  has  sections  with  significant  rate- 
limiting  properties  and  it  is  well-instrumented  for  comparison  with  computational  results.  The 
previous  attempt  to  model  the  Dover  site  was  a  two  dimensional  one  [2 1  ] ;  in  this  research  a  3D 
model  of  the  same  site  will  be  implemented. 

1.2  SPECIFIC  PROBLEM 

The  objective  of  this  work  is  to  mathematically  model,  in  three  dimensions,  aquifer  remedi¬ 
ation  by  continuous  and  pulsed  pumping  when  contaminant  transport  is  affected  by  macroscale 
sorption  and  diffusion.  This  research  will  extend  the  work  of  Herman  [21] 

1.3  RESEARCH  OBJECTIVES 

The  specific  objectives  of  this  research  are  to: 

1.  Install  and  test  the  Ground-water  modeling  System/FEMMATER  package  as  furnished  by 
the  Army  Corps  of  Engineers’  Waterways  Experiment  Station  on  AFIT  computational 
facilities. 

2.  In  a  three  step  process,  apply  GMS/FEMWATER  to  a  conceptual  site  developed  from  an 
actual  field  site  at  Dover  AFB,  Delaware: 

i  Conduct  a  moment  analysis  validation  of  the  numerical  formulation. 

ii  Determine,  through  numerical  simulation,  qualitative  aspects  of  remediation  by 
continuous  pumping  of  this  conceptual  site  with  special  attention  to  tailing  and  the 
effect  of  retardation  factor. 

iii  Determine  qualitiative  aspects  of  rebound  due  to  pulse  pumping. 


8 


2.  THE  GMS/FEMWATER  PACKAGE  -  AN  OVERVIEW 

This  section  will  explain  the  strengths  of  the  GMS  graphical  user  interface  and  the  FEMWA- 
TER  FORTRAN  program.  This  information  will  probably  be  dated  upon  publication  since  the 
software  in  question  is  constantly  being  revised  and  improved. 

2.1  FEATURES  OF  THE  GROUND- WATER  MODELING  SYSTEM  (GMS) 

The  Department  of  Defense  Ground- Water  Modeling  System  (GMS)  is  a  comprehensive 
graphical  user  environment  for  numerical  modeling.  It  was  developed  by  Brigham  Young  Uni¬ 
versity’s  Engineering  Computer  Graphics  Laboratory  in  cooperation  with  the  U.S.  Army  Corps 
of  Engineers  Waterways  Experiment  Station.  GMS  was  originally  intended  for  ground-water 
modeling  applications  and  currently  includes  specialized  interfaces  to  some  of  the  more  popu¬ 
lar  ground-water  modeling  programs  such  as  MODFLOW,  MT3D,  and  FEMWATER.  Flowever, 
GMS  is  written  in  a  general  fashion  so  that  it  can  be  used  as  a  platform  for  any  type  of  two  or 
three-dimensional  numerical  modeling. 

2.1.a  OPERATIONAL  MODULES 

The  GMS  screen  is  shown  in  figure  1 .  The  large  central  space  is  the  Graphics  Window. 
On  the  far  left  is  the  Tool  Palette-  a  collection  of  pushbutton  icons.  The  pull-down  Menu  Bar 
runs  across  the  top.  Finally,  the  Edit  Window  is  under  the  Graphics  Window  (to  the  right  of  the 
Tool  Palette). 

The  interface  for  GMS  is  divided  into  nine  separate  modules.  GMS  provides  a  module  for 
each  of  the  basic  datatypes  it  supports  such  as  Triangulated  Irregular  Network,  Borehole,  Solid, 
2D  and  3D  mesh,  or  3D  Scatter  Point.  As  the  user  switches  from  one  module  to  another,  the  Tool 
Palette  and  the  menus  change.  This  allows  the  user  to  focus  only  on  the  tools  and  commands 
relevant  to  the  data  type  the  user  will  manipulate  in  the  modeling  process.  The  user  can  switch 
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Figure  1.  The  GMS  Screen 
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instantaneously  from  one  module  to  another  to  facilitate  the  simultaneous  use  of  several  data 
types  when  necessary. 

2.1. b  DAIASETS 

The  interface  to  each  of  the  separate  modules  is  designed  in  a  consistent  fashion.  Once 
the  user  becomes  familiar  with  the  interface  to  one  of  the  modules,  the  other  modules  can  be 
used  immediately  with  little  further  training.  To  help  provide  a  consistent  interface,  GMS  uses 
the  concept  of  generic  data  sets.  A  data  set  is  a  set  of  scalar  or  vector  values  associated  with 
an  object.  Each  data  set  can  be  either  steady  state  or  transient  (multiple  values  representing  the 
data  values  at  different  points  in  time).  All  the  different  types  of  data  groups  have  associated 
lists  of  scalar  data  sets  and  vector  data  sets.  Each  set  has  a  single  vector  or  scalar  value  for  each 
node,  cell,  or  data  point. 

Data  sets  can  be  used  to  represent  a  variety  of  different  types  of  information.  They  can 
be  imported  from  a  file  or  they  can  be  created  by  interpolating  from  a  group  of  scattered  data 
points.  In  some  cases  it  is  necessary  to  perform  mathematical  operations  on  data  sets.  The 
user  can  accomplish  this  in  GMS  using  a  ’’data  set  calculator”.  For  example,  to  compare  the 
difference  in  the  solutions  from  two  separate  simulations  on  a  finite  difference  grid,  the  two 
solutions  can  be  input  as  data  and  the  data  calculator  used  to  compute  the  absolute  value  of  the 
difference  between  the  two  data  sets.  The  resulting  data  set  can  be  contoured  or  used  to  display 
iso-surfaces  just  like  any  other  data  set. 

2.1. C  VISUALIZATION 

GMS  provides  a  large  number  of  visualization  tools.  The  user  can  display  all  objects 
in  a  3D  oblique  view  and  rotate  them  interactively.  He  or  she  can  also  use  hidden  surface 
removal  and  color  shading  with  a  light  source  to  generate  highly  realistic  images.  In  addition, 
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the  modeler  can  use  contours  and  color  fringes  to  display  the  variation  of  input  data  or  computed 
data.  Starting  with  3D  meshes  and  grids,  he  or  she  can  generate  cross  sections  and  iso-surfaces. 
Given  an  object  with  an  associated  transient  data  set,  the  user  can  generate  animation  sequences 
showing  data  evolution. 

The  interface  to  the  visualization  tools  in  GMS  is  consistent  for  each  of  the  supported  data 
types.  The  dialogs  and  commands  used  for  visualization  are  identical  in  each  module. 

The  GMS  Reference  Manual  [9]  provides  complete  information  on  all  of  GMS’s  features. 

2.1.d  GRID  GENERATION 

Two  very  useful  grid-generation  features  in  GMS  are  the  adaptive  tessellation  2D  mesh 
construction  tool  and  the  3D  borehole-to-mesh  tool.  The  former  allows  the  user  to  generate 
2D  meshes  which  are  then  projected  onto  the  plan  geometry  of  the  site  in  question  to  give  a 
3D  mesh.  Adaptive  tessellation  is  advantageous  when  there  are  interior  regions  around  which 
finite  elements  must  be  wrapped.  When  a  boundary  polygon  encloses  these  interior  regions,  and 
adaptive  tessellation  is  used,  the  node  spacing  on  the  boundary  of  the  input  polygon  is  used  to 
determine  the  element  sizes  on  the  interior.  On  the  other  hand,  the  borehole-to-mesh  technique 
is  particularly  useful  when  the  site  stratification  character  is  provided  in  the  form  of  borehole 
data.  The  user  is  provided  the  means  to  connect  between  partitions  in  the  boreholes  and  directly 
fill  these  with  finite  elements,  making  for  rapid  visually-intuitive  mesh  generation. 

2.2  MATHEMATICAL  FORMULATION  BEHIND  FEMWATER 

Before  proceeding,  the  reader  should  note  that  the  essential  unknowns  in  the  unsteady 
equations  for  fluid  flow  and  contaminant  transport  are  hydraulic  head  h  and  contaminant  con¬ 
centration  C.  Transport  velocity  is  a  function  of  h.  As  will  be  detailed  below,  the  equation  gov¬ 
erning  fluid  flow  and  that  governing  contaminant  transport  are  coupled  through  an  unsteady 
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term  ^  and  through  the  dependence  of  water  density  and  dynamic  viscosity  on  contaminant 
concentration.  If  the  density  and  dynamic  viscosity  of  the  water  are  independent  of  contaminant 
concentration,  and  if  the  problem  under  consideration  is  steady  =  0),  the  equations  become 
uncoupled.  Then,  the  velocity  field  can  be  computed  independently  of  the  concentration  dis¬ 
tribution.  Once  this  velocity  field  is  determined,  it  enters  the  transport  equation  as  a  spatially- 
variable  coefficient.  The  transport  equation  can  then  be  solved  based  on  this  pre-computed  ve¬ 
locity  field.  The  following  two  sections  describe  the  governing  equations  in  FEMWATER  for 
flow  and  transport. 

2,2.a  GOVERNING  EQUATIONS  FOR  FLOW 

The  governing  equations  for  flow  through  aquifer  material  are  based  on  continuity  of  fluid, 
continuity  of  solid,  consolidation  of  the  media,  and  the  equation  of  state  [38,  pg  2] .  As  used  in 
FEMWATER,  these  governing  equations  for  flow  are: 


—F^  =  V-  [k-  (v/i  +  -^V2^1  +^q 


Pq  dt  L  V  Po  J\  Po 

where  p  is  the  fluid  density,  Pq  is  the  reference  water  density,  F  is  the  storage  coefficient, 
h  is  the  reference  pressure  head,  and  t  is  time.  Here  the  hydraulic  conductivity  tensor  K 
with  p  being  the  dynamic  viscosity  of  the  fluid,  g  the  acceleration  due  to  gravity  and  k  the 
intrinsic  permeability  tensor.  2  is  the  elevation  head,  /?*  the  density  of  the  injected  fluid,  and  q 


the  internal  source/sink. 


Once  the  flow  equation  is  solved  for  the  h-field,  the  Darcy  velocity  vector  for  density- 
dependent  flow  is: 

v  =  -K-  + 

The  density  and  dynamic  viscosity  are  functions  of  contaminant  concentration  and  are 


assumed  to  take  the  following  form: 


f  =ai+a2C  +  0302  + 04(^3 
—  ft5  +  (10(7  4-  (17(72  +  ClgC^ 

where  /.iq  is  the  reference  water  dynamic  viscosity,  C  is  the  contaminant  concentration 
(M/L3)  and  ai,  02, . . 07,  ag  (L^/M)  are  the  parameters  that  are  used  to  describe  the  concentra¬ 
tion  dependence  of  water  density  and  dynamic  viscosity.  Clearly,  if  02  =  03  =  04  =  ae  =  07  = 
ag  =  0  then  the  flow  equation  will  be  independent  of  contaminant  concentration  and  therefore 
uncoupled  from  the  transport  equation  (8). 

The  initial  conditions  for  the  flow  equations  are  stated  as: 

h  =  hi  (x,  z)  in  R 

where  R  is  the  region  of  interest  and  hi  is  the  prescribed  initial  condition  for  hydraulic 

head. 

In  order  to  describe  boundary  conditions,  we  must  first  define  some  nomenclature; 

n  is  the  outward  unit  vector  normal  to  the  boundary;  {xb,yi,,  Zb)is  the  spatial  coordinate  on 
the  boundary;  hd,  qn,  Qc  are  the  Dirichlet  functional  value  for  hydraulic  head,  Neumann  flux,  and 
Cauchy  flux,  respectively;  Bn,  Be,  By  are  the  Dirichlet,  Neumann,  Cauchy,  and  variable 
boundaries,  respectively.  Then  the  boundary  is  the  union  of  Bd,  Bn,  Be,  and  By.  Next,  hp  and 
qp  and  are  the  allowed  ponding  depth  and  the  throughfall  of  precipitation,  respectively,  on  the 
variable  boundary;  hm  is  the  allowed  minimum  pressure  on  the  variable  boundary;  and  qe  is  the 
allowed  maximum  evaporation  rate  (=  potential  evaporation)  on  the  variable  boundary;  for  a 
boundary  on  a  river  (a  common  boundary  for  flow  problems),  Kr  is  the  hydraulic  conductivity 
on  the  river  bottom  sediment  layer,  ha  is  the  thickness  of  the  river  bottom  sediment  layer,  and 
hfi  is  the  depth  of  the  river  bottom  measured  from  the  river  surface. 

Five  types  of  boundary  conditions  can  be  specified  for  the  flow  equations  depending  on 
the  physical  location  of  the  boundaries.  These  boundary  conditions  are  stated  as: 
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Dirichlet  Boundary  Conditions: 


h  =  hd  {xb,yb,Zb,t) 


on  Bd 

Neumann  Boundary  Conditions: 


qn{xb,yb,Zb,t) 


on  B„ 


Cauchy  Boundary  Conditions: 


-n  K 


^Vh  +  Vz 
,  P 


<ic{xb,yb,zb,t) 


on  Be 

Variable  Boundary  Conditions  -  During  Precipitation  Period: 


h  =  hp{xb,yb,Zb,t) 


on  Bt, 


or 


on  B^ 


-n  K 


qpixb,yb,zb,t) 


Variable  Boundary  Conditions  -  During  Non-Precipitation  Period: 


3 


on  B,, 


or 


h  —  hp  {xb,ybt  Zbi i) 


4 
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h  =  h,n  {Xb,yb,ZbJ) 


5 


on  By 


or 


on  B^, 


-nK 


Qe  {Xb-,yb,Zb,t) 


6 


River  Boundary  Conditions: 


on  Br 


-n  K 


Po 


Vh  +  Vz 


1 


Note  that  only  one  of  Eqs.  (2)  through  (6)  is  utilized  at  any  point  on  the  variable  boundary 


at  any  time. 
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2.2.b  GOVERNING  EQUATIONS  FOR  TRANSPORT 

The  governing  equations  for  material  transport  through  ground-water  systems  are  derived 
based  on  the  laws  of  continuity  of  mass  and  flux.  The  major  processes  that  are  included  are  ad- 
vection,  dispersion/diffusion,  decay,  adsorption,  biodegradation  through  both  liquid  and  solid 
phases,  the  compressibility  of  media,  as  well  as  sources  and  sinks.  Letting  C  be  the  dis¬ 
solved  concentration  and  S  be  the  adsorbed  concentration,  the  governing  equation  of  the  spatial- 
temporal  distribution  of  dissolved  concentrations  can  be  stated  as  follows; 

^  '"'’If  ^  ^  ^ 


+  {9C  +  p,S)  -  {OK^C  +  Pt,KsS)  + 


m 


at  p 


30 
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p-  \  at  p  \poj  dt 
where  6  is  the  dimensionless  moisture  content,  is  the  bulk  density  of  the  medium,  V  is 

the  darcy  velocity,  D  is  the  matrix  of  dispersion  coefficients,  A  is  the  decay  constant  (l/T),  Kyj 

is  the  first  order  biodegradation  constant  through  dissoved  phase  (l/T),  Ks  is  the  first  order 

biodegradation  constant  through  adsorbed  phase  (l/T),  q  is  the  internal  source/sink,  and  m  is 

external  source/sink  rate  per  medium  volume  ((M/L^)/T).  Equation  (8)  above  involves  two 

unknowns  C  and  S,  so  for  closure,  constitutive  relationships  must  be  posed.  In  FEMW\TER, 

the  following  empirical  relationships  are  used: 

For  the  linear  isotherm: 


S^KaC 

For  the  Langmuir  isotherm: 

^  “  l+KC 

For  the  Freundlich  isotherm: 
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S  =  KC^ 


where  Ka  is  the  distribution  coefficient  (L^/M),  Smaxis  the  maximum  concentration  per¬ 
mitted  in  the  medium  in  the  Langmuir  nonlinear  isotherm,  K  is  the  coefficient  in  the  Langmuir 
or  Freundlich  nonlinear  isotherm,  and  n  is  the  power  index  in  the  Freundlich  nonlinear  isotherm. 
Note  that  all  three  of  these  models  relate  dissolved  concentration  to  adsorbed  concentration  by 
using  the  local  equilibrium  assumption. 

Assuming  that  initially  the  dissolved  concentrations  are  known  throughout  the  region  of 
interest: 

C  =  Ci  (x,  z)  in  R 

where  Ci  is  the  initial  concentration  and  R  is  the  region  of  interest.  Initial  concentrations 
for  the  dissolved  concentrations  may  be  obtained  from  field  measurements. 

The  specification  of  boundary  conditions  is  a  difficult  and  intricate  task  in  transport  mod¬ 
eling.  From  the  dynamic  point  of  view,  a  boundary  segment  may  be  classified  as  either  flow¬ 
through  or  impervious.  From  a  physical  point  of  view,  it  is  a  soil-air  interface,  or  soil-soil 
interface,  or  soil-water  interface.  From  the  mathematical  point  of  view,  it  may  be  treated  as  a 
Dirichlet  boundary  on  which  the  total  analytical  concentration  is  prescribed,  a  Neumann  bound¬ 
ary  on  which  the  flux  due  to  the  gradient  of  total  analytical  concentration  is  known,  or  a  Cauchy 
boundary  on  which  the  total  flux  is  given. 

An  even  more  difficult  mathematical  boundary  is  the  variable  boundary  conditions  on 
which  the  boundary  conditions  are  not  known  a  priori  but  are  themselves  the  solution  to  be 
sought.  In  other  words,  on  the  mathematically  variable  boundary,  either  Neumann  or  Cauchy 
conditions  may  prevail  and  change  with  time.  Which  condition  prevails  at  a  particular  time  can 
be  determined  only  in  the  cyclic  processes  of  solving  the  governing  equations,  iterating  time- 
step  by  time-step  between  solution  fields  and  boundary  conditions[38,  pg  95]  . 
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Whatever  point  of  view  is  chosen,  all  boundary  conditions  eventually  must  be  transformed 
into  mathematical  equations  for  quantitative  solutions.  Thus,  we  will  specifj'  the  boundary  con¬ 
ditions  from  the  mathematical  point  of  view  in  concert  with  dynamic  and  physical  considera¬ 
tions.  The  boundary  conditions  imposed  on  any  segment  of  the  boundary  are  taken  to  be  ei¬ 
ther  Dirichlet,  Neumann,  Cauchy,  or  variable.  Thus,  the  global  boundary  may  be  split  into  four 
parts,  Bd,Bn,Bc,  and  By  denoting  Dirichlet,  Neumann,  Cauchy,  and  variable  boundaries,  re¬ 
spectively.  The  conditions  imposed  on  the  first  three  types  of  boundaries  are  given  as: 
Prescribed  Concentration  (Dirichlet)  Boundary  Conditions: 

C  =  Cd  {xb,  Vb,  Zb,  t)on  Bd 

Neumann  Boundary  Conditions: 

n  •  (-^D  •  VC)  =  qn  {xb,yb,Zb,t)on  Bn 

Cauchy  Boundary  Conditions: 

n  •  (VC  -  6>D  •  VC)  =  qc{xb,yb,Zb,t)on  Be 

where  Cd  is  the  prescribed  concentration  on  the  Dirichlet  boundary  Bd.  {xb,yb,  Zb)is  the 
spatial  coordinate  on  the  boundary,  n  is  an  outward  unit  vector  normal  to  the  boundary,  6  is  the 
moisture  content,  qn  is  the  prescribed  gradient  flux  through  the  Neumann  boundary  J5„,  and  qc 
is  the  prescribed  total  flux  through  the  Cauchy  boundary  Be. 

The  conditions  imposed  on  the  variable-type  boundary,  which  is  normally  the  soil-air  in¬ 
terface  or  soil-water  interface,  are  either  the  Neumann  with  zero  gradient  flux  or  the  Cauchy 
with  given  total  flux.  The  former  is  specified  when  the  water  flow  is  directed  out  of  the  region 
from  the  far  away  boundary,  whereas  the  latter  is  specified  when  the  water  flow  is  directed  into 
the  region.  This  type  of  variable  condition  would  normally  occur  at  flow-through  boundaries. 
Written  mathematically,  the  variable  boundary  condition  is  given  by 
n  •  (VC  -  6>D  •  VC)  =  n  •  VCy  {xb,yb,Zb,t)on  5^,  if  n  ■  V  <0 
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n  •  (-0D  •  VC)  =  0  on  if  n  •  V  >0 

where  V  is  the  Darcy  velocity,  C„  is  the  specified  concentration  of  water  through  the 


variable  boundary  and  is  the  variable  boundary. 

An  issue  to  be  aware  of  is  the  possibility  of  ill-posing  the  transport  problem  which  may 
result  in  a  non-unique  solution.  For  a  simplified  mathematical  explanation  see  appendix  A. 

2.2.C  NUMERICAL  FORMULATION 

The  interested  reader  may  see  [38,  pg  97-123]  for  details  on  the  numerical  formulation  of 
FEMWATER.  However,  some  of  the  notable  numerical  features  of  the  code  are: 

1.  Spatial  discretization  is  of  the  finite  element  variety.  Time  discretization  is  by  finite 
differences. 

2.  The  flow  module  in  FEMWATER  has  the  following  features: 

A.  Galerkin  finite  element  method  for  spatial  discretization 

B.  To  linearize  the  matrix  equation,  the  Picard  method  is  used  instead  of  the  Newton- 
Raphson  which  would  result  in  an  asymmetric  matrix. 

C.  In  solving  the  linearized  matrix  equations,  direct  methods  aren’t  practical  in  dealing 
with  large  3-dimensional  problems.  Instead,  three  iteration  modules  are  available  to 
the  user: 

*  Successive  point  iteration 

*  Polynomial  preconditioned  conjugate  gradient 

*  Incomplete  Cholesky  preconditioned  conjugate  gradient 

D.  To  handle  the  mass  matrix  resulting  from  the  storage  term,  2  options  are  available 
(lumping  &  consistent) 

E.  In  approximating  the  time  derivatives,  2  options  are  available  (time-weighted 
difference  and  mid-difference). 

3.  The  transport  module  in  FEMWATER  provides  the  equation-solving  options  of  the  flow 
model  plus  three  finite  element  options: 

-  Galerkin 

-  Upstream  weighting 

-  Hybrid  Lagrangian-Eulerian 
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3.  CHARACTERISTICS  OF  THE  DOVER  AFB  SITE 


This  numerical  study  parallels  an  experimental  effort  conducted  by  Ball  et  al  [27]  of  the 
Johns  Hopkins  University  in  collaboration  with  the  University  of  Waterloo  and  the  Air  Force 
Research  Laboratory  (Tyndall  AFB,  Florida).  The  primary  objective  of  the  experiment  was  to 
determine  whether  there  are  discernible  advantages  to  a  pulsed-pumping  strategy  compared  to 
continuous  pumping  at  a  ”real-world”  site  where  the  ground  water  and  subsurface  solids  had 
long  been  contaminated  by  volatile  organic  contaminants  (VOCs).  The  underlying  scientific 
goal  was  to  identify  the  most  significant  sources  of  mass  transfer  rate  limitation  at  the  field  site 
and  to  develop  appropriate  conceptual  and  computational  models  to  describe  and  predict  the 
resulting  effects  on  aquifer  decontamination.The  experimentalists  chose  the  specific  site  for  its 
following  attributes; 

1.  Well-studied,  mildly  heterogeneous,  moderately  sorbing  sand/gravel  aquifer  impacted  for 
at  least  a  decade  with  chlorinated  volatile  organic  contaminants  (VOCs)  at  concentrations 
several  orders  of  magnitude  above  typical  practical  quantification  levels. 

2.  The  water  table  is  within  20  feet  of  ground  surface  (to  allow  sampling  by  pumps  at  the 
surface). 

3.  A  competent  (confining)  clay  aquitard  within  50  feet  of  ground  surface  to  form  the  bottom 
of  the  isolated  portions  of  the  aquifer. 

4.  Available  power  and  other  utilities,  ease  of  access,  and  security 

3.1  SITE  DESCRIPTION 

Figure  2  is  a  schematic  top  view  of  the  experimental  site.  Two  cells,  of  roughly  identical 
physical  dimension,  lie  beside  each  other.  By  design,  the  experimentalists  had  two  practically 
identical  cells  in  which  to  compare  the  remediative  performances  of  continuous  versus  pulsed 
pumping  -  each  technique  exclusively  applied  to  a  single  cell. 

The  cells  were  sealed  by  driving  sheet  piling  from  the  ground  surface  through  the  aquifer 
layers  into  the  underling  clay  aquitard.  Each  cell  was  then  cored  (borehole  sampled)  and  instru- 
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merited  with  various  types  of  wells.  Three  injection  and  three  extraction  wells  were  installed  at 
the  opposing  ends  of  both  cells.  These  were  fully  screened  and  penetrated  to  just  above  the  clay 
layer.  For  complete  details  on  the  experimental  procedure,  see  [27,  pg.  37-42] 
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Figure  2.  Schematic  of  Dover  AFB  Field  Site  (not  to  scale) 
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3.2  DEFINING  THE  MODEL  FOR  FEMWATER 


3.2.a  CONCEPTUAL  CELL 

Since  the  two  actual  cells  are  nearly  identical,  it  suffices,  for  modeling  purposes,  to  define 
a  single  "nominal”  cell  from  a  simple  average  of  the  two  cells’  dimensions.  This  conceptual 
cell  [21,  pg  9]  is  3.8  meters  wide  by  10  meters  long  by  13  meters  deep.  The  cell,  like  its  parents, 
is  heterogeneous  but  well  stratified  vertically.  A  simplified  conceptual  model  assuming  perfect 
stratification  was  developed  from  early  field  data.  Figure  3,  based  on  experimental  borehole 
data,  is  a  vertical  cross  section  schematic  of  the  soil  layers. 

Note  that  for  numerical  stability  reasons,  a  thin  layer  (0.15  m  thick)  of  fictitious  material, 
called  "Sandy  Loam”  was  inserted  between  the  Orange  Coarse  Sand  layer  and  the  Orange  Silty 
Clay  layer.  This  material’s  hydraulic  conductivity  was  halfway  between  those  of  the  two  adja¬ 
cent  materials.  This  addition  is  in  accordance  with  guideline  2  provided  in  section  3.2. b  below. 
Tables  1  and  2  show  each  material’s  physical  properties  as  used  in  the  numerical  simulation. 

The  four  Rhombus-shaped  symbols  below  the  rightmost  tip  of  the  concentration  profile  in 
figure  3  represent  the  vertical  positions  of  the  four  monitoring  points  in  the  geometric  center  (x  = 
5  m  and  y  =  1 .9  m)  of  the  test  cell.  The  placement  of  these  monitoring  points  was  in  accordance 
with  the  following: 

1 .  MPl  was  at  the  edge  of  the  clay  layer 

2.  MP2  corresponds  to  Herman’s  [21]  "central  monitoring  point”,  6  inches  above  the  clay 
layer 

3.  MP4  was  at  the  peak  of  preliminary  profiles  (supplied  by  [27]  )  of  contaminant 
concentration  versus  depth 

4.  MP3  was  approximately  half-way  between  MP2  and  MP4 

The  exact  positions  of  these  monitoring  points,  MPl  through  MP4,  are  shown  in  Table  3 
below: 
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SOIL  STRATA  &  MONITOR  POINTS 


Figure  3.  Soil  Layers  and  4  Monitoring  Points 


Table  1 .  Material  Properties  for  Sand  Layers 


Tan 

Medium 

Sand 

Gravel 

Tan 

Sand 

Gray 

Medium 

Sand 

Orange 

Medium 

Sand 

Orange 

Coarse 

Sand 

Distribution  Coefficient  (^) 

10-'* 

10-4 

10-4 

10-4 

10-4 

Bulk  Density  (^) 

1700 

1700 

1700 

1700 

1700 

Longitudinal  Dispersivity  (m) 

0.23 

0.23 

0.23 

0.23 

0.23 

Lateral  Dispersivity  (m) 

0.023 

0.023 

0.023 

0.023 

0.023 

Molecular  Diffusivity  (^) 

2.98-10-® 

2.98-10-® 

2.9810-® 

2.98- 10-® 

2.9810-® 

Tortuosity 

0.6 

0.6 

0.6 

0.6 

0.6 

Decay  Coefficient  (^) 

0 

0 

0 

0 

0 

Freundlich  N 

1 

1 

1 

1 

1 

Conductivity  (^) 

1.08-10-^ 

4.32-10-^ 

1.08-10-^ 

5.4-10-^ 

8.28-10-'-^ 

Moisture  Content 

.36 

.36 

.36 

.36 

.36 

Relative  Conductivity  (^) 

1 

1 

1 

1 

1 

Water  Capacity  (;^) 

0 

0 

0 

0 

0 

Table  2.  Material  Properties  for  Clay  Layers 


Orange 

Silty 

Clay 

Loam 

Black 

Silty 

Loam 

Sandy 

Loam 

(Fictitious) 

Distribution  Coefficient  (^) 

3.6-10-4 

1.94-10-2 

2.3-10-4 

Bulk  Density  (^) 

1300 

1300 

1500 

Longitudinal  Dispersivity  (m) 

0.23 

0.23 

0.23 

Lateral  Dispersivity  (m) 

0.023 

0.023 

0.023 

Molecular  Diffusivity  (^) 

2.98-10-® 

2.98-10-® 

2.98-10"® 

Tortuosity 

0.6 

0.6 

0.6 

Decay  Coefficient  (^) 

0 

0 

0 

Freundlich  N 

1 

1 

1 

Conductivity  (^) 

3.6-10-'^ 

3.6-10-^ 

ir^^ 

Moisture  Content 

.36 

.36 

.36 

Relative  Conductivity  (^) 

1 

1 

1 

Water  Capacity  (;^) 

0 

0 

0 

Table  3.  \fertical  Positioning  of  Central  Monitoring  Points 


Monitoring  Point 

MPl 

MP2 

MP3 

MP4 

Mesh  Node  Number 

1 

9334 

12445 

14519 

\fertical  Position 

-10m 

-9.8476m 

-9.4327m 

-8.8788m 
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3.2.b  FINITE  ELEMENT  MESH 

Figure  4  shows  the  three-dimensional  finite  element  mesh  used  for  this  simulation.  It 
consists  of 33 1 84  nodes.  The  corresponding  finite  elements  are  right  prisms  of  triangular  cross- 
section  (figure  5). 

Salient  features  of  the  mesh  are: 

1.  Fine  lateral  clustering  of  nodes  in  the  vicinities  of  the  three  injection  and  three  extraction 
wells  (figure  6). 

2.  J^rtical  clustering  in  accordance  with  the  two  FEMWATER  guidelines  [38,  pg  16]  calling 

A.  No  more  than  50%  increase  or  decrease  in  widths  of  adjacent  elements. 

B.  A  minimum  of  three  layers  of  elements  vertically  for  each  distinct  stratigraphic  unit, 
particularly  if  large  variations  of  hydraulic  conductivity  occur  in  adjacent  layers. 

C.  No  more  than  three  orders  of  magnitude  difference  in  hydraulic  conductivity  between 
adjacent  mesh  layers. 
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Figure  5.  Plan  View  of  3D  mesh 


Figure  6.  Mesh  Detail  in  the  Vicinity  of  Well 
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3.2.C  BOUNDARY  CONDITIONS 


In  imposing  boundary  conditions,  two  assumptions  were  made:  First,  that  the  lower  bound¬ 
ing  surface  of  the  cell  was  in  fact  a  confining,  impermeable  one.  Second,  that  there  was  no 
ponding  and  that  the  water  table  coincident  with  the  upper  surface  of  the  cell.  Then,  the  bound¬ 
ary  conditions  were  simply  zero  flux  of  flow  or  contaminant  on  all  six  sides  of  the  test  cell.  The 

3 

total  volumetric  pump  rate  for  three  injection  wells  was  Qt  =  0.192 

3.2.d  DISCRETIZATION  OF  WELLS 

Here,  the  challenge  was  to  model  an  injection  or  extraction  well  -  essentially  a  cylindri¬ 
cal  sieve  -  in  a  discrete  manner  consistent  with  our  finite  element  mesh.  As  shown  in  figure  7, 
Each  injection  and  extraction  well  was  modeled  by  a  series  of  sources  (or  sinks,  as  appropriate). 
It  would  have  been  an  over-simplification  to  simply  divide  the  total  injection  (or  extraction) 
flow  rate  Qt  by  the  total  number  of  well  nodes  at  the  injection  (or  extraction)  end  to  deter¬ 
mine  the  individual  source  (or  sink)  strength.  Had  this  been  done,  the  result  would  have  been 
that  each  source  or  sink  would  have  had  identical  strength  and  this  would  have  been  inconsis¬ 
tent  with  the  fact  that  each  source’s  strength  must  be  proportional  to  the  hydraulic  conductivity 
of  it’s  confining  material.  Since  each  source  was  sandwiched  between  materials  of  different 
conductivities,  it’s  respective  strength  had  in  turn  to  be  different.  Referring  to  figure  8,  there 
were  several  underlying  assumptions  that  had  to  be  made.  First,  that  the  water  acting  under 
the  pressure  gradient  ^  induced  by  the  pumping  process,  moved  from  the  injection  wells  to 
the  extraction  wells  in  parallel  horizontal  planes  albeit  at  different  velocities.  Second,  that  this 
pressure  gradient  ^  was  then  constant  across  stream  layers.  These  assumptions  were  consis¬ 
tent  with  observations  made  in  seminal  works  such  as  [2] .  Then,  the  Darcy  equation  [11,  pg  16] 
Qt  —  ^  H  ^  held  not  only  for  the  whole  cell  (of  width  W  and  height  H),  but  also  for 
each  individual  stream  layer  (of  width  w  and  height  6s j).  The  stream  layer  width  w,  depended 
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on  the  relative  strengths  of  adjacent  wells.  Referring  to  figure  7,  The  lateral  hydraulic  zone  of 
influence  of  each  well  was  in  direct  proportion  to  that  well’s  volumetric  flow  strength.  For  sim¬ 
plicity,  each  well  was  assigned  equal  volumetric  flow  strength.  Then  Qj  =  Kj  ■  w  ■  6sj  •  ^ 
could  be  written  for  each  layer,  where  6sj  was  the  thickness  of  the  stream  layer  in  question  and 
kj  a  weighted  hydraulic  conductivity  relevant  to  6sj  since  6s j  spanned  adjacent  finite  element 
mesh  layers,  of  differing  hydraulic  conductivity.  If  the  total  volumetric  flow  rate  of  an  injec¬ 
tion  well  was  Qt,  then  the  volumetric  flow  rate  through  stream  layer  6s j  was  simply  Qj  where 
Qt  =  E  =  E  ^  ^  ■  S  S  ^  ^  summation 

being  over  all  stream  layers  j.  Note  that  in  what  follows,  i  indexes  mesh  layer  and  j  indexes 
stream  layer.  As  shown  in  figure  8,  the  stream  layers  overlapped  mesh  layers  A  vertical  series  of 
point  flow  sources  of  volumetric  strengths  Qj  was  placed  along  the  axis  of  the  well  in  question, 
the  sum  of  strengths  Qt  producing  the  overall  volumetric  flow  rate  of  the  well.  It  then  remained 
to  determine  the  relative  weighting  of  each  of  these  point  sources,  consistent  with  the  local  hy- 
draulic conductivities ifj.  Well,(5j  =  Kj-w-6sj-§  butru-f  = 

so  Qj  =  •  Qt-  Now  a  third  assumption  had  to  be  made,  concerning  our  choice  of 

the  6si,  the  stream  layer  thicknesses.  It  was  assumed  that  within  a  given  layer  of  mesh  elements, 
between  two  adjacent  sources,  the  sources  ’’shared”  the  layer  of  elements  equally.  Specifically, 
each  source’s  hydraulic  zone  of  influence  extended  to  the  centerline  of  the  adjacent  layer  of  el¬ 
ements.  In  the  exceptional  case  that  a  given  layer  of  mesh  elements  was  bounded  on  only  one 
edge  by  a  source,  then  that  source’s  hydraulic  zone  of  influence  encompassed  the  entire  width 
of  the  semi-bounded  element  layer.  This  is  depicted  schematically  in  figure  8.  Each  6si  was  the 
sum  of  the  local  6zti  and  6zbi.  It  was  assumed  6zti  and  6zbi-.i  were  equal  within  each  mesh 
layer  6zi.  In  other  words  6zti  =  6zbi-i  =  For  which,  Qi  =  ^ j  •  Qt 
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Figure  8.  Source  Placement  and  Streamsurface  Geometry 
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3.2.e  CURVES  OF  SOIL  HYDRAULIC  PROPERTIES 


The  FEMWATER  model  calls  for  the  following  three  curves  of  soil  hydraulic  properties 
versus  pressure  head  h\ 

1 .  Moisture  content  {6w) 

2.  Water  capacity  (^) 

3.  Relative  conductivity  (isTr) 

Since  the  test  cell  was  taken  to  be  saturated,  these  three  properties  are  constant  with  h 
[11,  pg  39]  .  However,  according  to  Zakikhani  [41]  ,  it  is  prudent  to  include  portions  of  the 
properties’  curves  for  unsaturated  regions  of  h  since,  during  the  transient  course  of  computing 
a  saturated  solution,  it  is  possible  for  the  computed  solution  to  dip  into  the  unsaturated  domain. 


Not  properly  representing  this  domain  through  the  soil  properties  curves  may  lead  to  instability. 
Hence,  computation  of  curves  was  done  using  the  full  Van  Genuchten  models  [36] . 

First,  the  so-called  ’’effective  moisture  content”,  6^  must  be  defined  —  in  which  a,  13,  and 
■y  =  1  —  i  are  constant  properties  of  the  material  in  question  : 

g(h)-l  +  if  h<0 

\  1  i/  /i>0 

Then,  with  9s  and  9r  (also  material  constants)  defined  as  ’’saturation  moisture  content”  and 
’’residual  moisture  content”  respectively,  moisture  content  is  defined  as: 


9-u}  —  Gj-  9e{9s  —  9r) 

Next,  ^{h)  is  given  by 

f  -  [1  +  ^  if  h<0 

dh^  ’  \  0  if  h>0 

Then, 

And  the  Van  Genuchten  model  for  relative  conductivity  (Kr)  is  given  by: 
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Kr{h)  =  ^/^h)  [l  -  (l  -  Oe{h)^) 


Using  data  provided  in  Van  Genuchten’s  original  paper  [36]  ,  as  well  as  sample  data  provided 
in  the  FEMW\TER  manual  [38]  as  a  basis,  values  for  the  the  parameters  a,  (3,  9s,  and  9r  were 
approximated.  These  are  listed  in  Table  4  below. 

Choosing  sand  for  representative  purposes,  figures  9,  10,  and  11  show  9e,  and  Kr 
respectively  versus  pressure  head  h.  Note  that  1000  points  were  used  for  each  curve  to  provide 
sufficient  resolution  in  the  vicinity  of  the  inflections  at  /i  =  0. 


Table  4.  Van  Genuchten  Parameters  for  Soil  Hydraulic  Properties 


9s 

9r 

a(;^) 

/? 

SAND 

0.43 

0.045 

14.5 

2.68 

SILTY  CLAY  LOAM 

0.43 

0.089 

1 

1.23 

SILTY  LOAM 

0.45 

0.067 

2 

1.41 
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Figure  10.  Water  Capacity  vs.  Pressure  Head  (Sand) 
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Figure  11.  Relative  Conductivity  vs.  Pressure  Head  (Sand) 
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4.  NUMERICAL  SIMULATIONS 


In  accordance  with  the  research  objectives  stated  in  section  1.3;  namely,  to  validate  the 
numerical  formulation  and  determine  cjualitative  aspects  of  pump-and-treat  remediation,  the 
following  simulations  were  run: 

1.  The  flow  solution  showing  the  velocity  field  as  vectors.  In  the  subsequent  computations, 
this  computed  velocity  field  formed  the  basis  for  all  transport  simulations.  As  explained 
in  section  2.2  above,  if  the  flow  is  steady-state,  the  velocity  field  enters  into  the  transport 
equation  as  a  spatially- variable  coefficient.  With  the  velocity  field  in  hand,  the  transport 
equation  can  then  be  solved. 

2.  The  transport  solution  resulting  from  the  injection  of  a  unit  concentration  of  contaminant 
for  200  hours.  This  solution  allowed  a  moment  analysis  cross-check. 

3.  Another  transport  solution  for  a  200  hour  injection  pulse,  as  above,  but  with  retardation 
factor  equal  to  one  for  all  soil  layers.  This  was  to  examine  the  effect  of  retardation  on 
contaminant  transport. 

4.  The  transport  solution  resulting  from  the  continuous  pumping  of  the  test  cell  with  clean 
water,  starting  with  a  prescribed  initial  contaminant  distribution.  This  was  to  spatially 
examine  how  the  aquifer/aquitard  composite  was  releasing  or  retaining  contaminant. 

5.  The  transport  solution  resulting  from  continuous  pumping  (as  above)  for  3600  hours, 
followed  by  pump  stand-down  for  an  additional  1900  hours,  then  resumption  of  pumping. 
This  was  to  examine  the  effect  of  diffusive  contaminant  rebound  from  the  aquitard  sublayer. 

6.  The  transport  solution  resulting  from  the  continuous  pumping  of  the  test  cell  with  clean 
water,  starting  with  an  initial  contaminant  distribution  provided  by  the  field  experimenters 
[27]  .  This  was  to  determine  how  well  the  FEMWATER  simulation  mimicked  data  from 
the  Dover  AFB  remediation  experiment. 

7.  The  transport  solution  resulting  from  the  continuous  pumping  of  the  test  cell  with  clean 
water,  starting  with  a  prescribed  initial  contaminant  distribution  provided  by  Herman  [21, 
pg  12]  .  This  was  to  compare  with  another  numerical  simulation  of  the  Dover  AFB 
remediation  experiment. 

Sections  4.1  through  4.7  detail  the  seven  simulations  listed  above. 
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4.1  COMPUTATION  OF  FLOW  SOLUTION 


A  sample  input  file  is  given  in  appendix  B.l.  The  computed  velocity  field  is  shown  in 
figures  12, 13,  and  14.  For  clarity,  only  every  third  computed  velocity  vector  is  shown  in  plot  12. 
Note  that  in  figure  14  the  large  arrows  depicting  vector  velocity  in  the  immediate  vicinity  of  the 
sink  appear  to  be  directed  away  from  the  sink.This  is  an  anomaly  of  the  plotting  software.  They 
should  be  oriented  so  that  there  heads  converge,  or  ’’touch”  at  the  sink.  However,  this  would 
obscure  detail  in  that  area,  hence  their  "backwards”  orientation.  There  are  several  observations 
that  can  be  made  immediately; 

1.  The  graphical  disparity  between  flow  admittance  in  the  various  aquifer  and  aquitard 
materials. 

2.  The  fact  that  the  fluid  ~  albeit  at  different  velocities  -  moves  in  parallel  layers,  supporting 
shear-generated  mixing. 

3.  Figure  14  is  a  close  up  of  the  flow  field  in  the  vicinity  of  a  sink  node.  It’s  remarkable  (but 
not  really  surprising!)  that  the  flow  field  is  very  three-dimensional  in  this  area  but  very 
ordered  and  laminar  near  the  center  of  the  test  cell  (figure  13). 

Finally,  figure  1 5  shows  the  vertical  distribution  of  longitudinal  velocity  at  the  center  of  the 
cell.  Only  every  third  computed  velocity  vector  is  shown  in  the  plot.  Once  again,  the  division 
of  flow  between  aquifer  and  aquitard  layers  is  very  evident. 
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Figure  12.  \felocity\fector  Field 
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TAN  MEDIUM  SAND  K  =  3  x  lO?  m/s 
GRAVEL  /  TAN  SAND  K  =  1 .2  x  10-^  m/s 

GRAY  MEDIUM  SAND  K==  3  x  10‘7  m/s 


ORANGE  MEDIUM  SAND 
K=  1.5  X  10-®  m/s 


ORANGE  MEDIUM  /  COARSE  SAND 
K=2.3  X  10-®  m/s 

ORANGE  SILTY  CLAY  LOAM  K=  10  ^  m/^ 
E3LACK  SILTY  LOAM  K  =  10'S  m/s 


Figure  15.  \fertical  Distribution  of  Longitudinal  \felocity  (m  /  hr) 
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Normally,  the  validity  of  the  flow  solution  would  be  aeertained  by  comparing  the  distri¬ 
bution  of  computed  pressure  head  (from  which  the  velocity  field  is  derived)  with  the  pressure 
head  distribution  measured  in  the  field.  However,  the  latter  was  unavailable  at  the  time  of  writ¬ 
ing.  Nevertheless,  the  qualitative  behaviour  indicated  in  the  observations  above  was  consistent 
with  expectations.  A  moment  analysis  (following  section)  of  the  results  provided  additional 
validation. 
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4.2  MOMENT  ANALYSIS  CROSS-CHECK 

The  cross-check  was  an  order-of-magnitude  comparison  between  two  calculations  of  mean 
residence  time  of  contaminant  in  transit  between  the  injection  wells  and  the  extraction  wells. 
The  first  (Tres)  is  based  solely  on  the  physical  characteristics  of  the  test  cell.  The  second  (T^es) 
is  based  on  a  moment  analysis  of  the  computed  curve  of  mass  recovery  rate  versus  time  at  the 
extraction  well.  The  first  calculation  is  as  follows: 

1  A  mean  velocity  is  calculated  at  the  center  of  the  cell  over  the  depth  of  the  cell: 


fm.ean  — 


£^1  ■  (I’M  -  hi) 


hhr  —  hi 


0.003882  m/hr 


where  Vi  and  are  computed  velocities. 


2  A  transport  length,  L,  is  given  by  the  longitudinal  distance  between  the  line  of  injection 
wells  and  the  line  of  extraction  wells: 

L  —  ^extvcLction  ^injection  9.6654  0.3866  9.2T88  m 


3  A  composite  retardation  factor  Ream  is  given  by  the  expression: 

^  Capacity  for  solute  in  all  zones 

Capacity  for  solute  in  mobile  zones 

n  Eili  porosityi  ■  depthi  •  Riayer,i  _  ^  o  o 

xXcOTTl  . '  ““  -1.  ^  ^ 

Z^i=i  porosityi  ■  depthi  ■  Riayer,i 

where  Riayer,i  is  the  retardation  coefficient  associated  with  layer  i  -  of  thickness  depthi. 

4  Finally,  the  mean  residence  time  based  on  mean  velocity  Kneon  and  composite  retardation 
factor  Room  is 


Tres  — 


-Rrnm  =  29161  hOUtS 


The  second  calculation,  is  based  on  the  moment  analysis  [14,  pg  1575-1585]  of  the 
contaminant  breakthrough  curve  at  the  extraction  well.  The  idea  is  to  compute,  using  the  sim¬ 
ple  principle  of  moments,  the  "center  of  gravity”  of  the  area  under  the  curve  of  mass  recovery 
rate  versus  time.  In  the  present  case  the  observation  location  is  the  line  of  extraction  wells  and 
is  a  measure  of  the  mean  residence  time  of  contaminant  in  the  test  cell.  The  actual  numeri¬ 
cal  simulation  proceeded  as  follows:  A  unit  concentration  of  contaminant  was  injected  into  the 
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test  cell  for  a  known  length  of  time  (200  hours).  This  was  the  test  pulse.  The  introduction  of 
contaminant  was  then  ceased  but  the  flushing  cycle  with  clean  water  (at  the  original  pumping 
rate)  was  maintained.  Throughout  this  process,  the  mass  collection  rate  at  the  extraction  well 
was  monitored.  This  is  defined  as  the  iimer  product  of  two  vectors;  namely,  contaminant  con¬ 
centrations  at  the  sink  nodes  and  corresponding  sink  flow  rates,  i.e.  Ci-Qi.  The  curves 

depicting  injection  and  extraction  mass  rates  respectively  are  shown  schematically  in  figure  16 
below.  The  rectangular,  injection  pulse  is  on  the  left.  The  bell-shaped  ’’breakthrough”  curve, 
describing  extraction  mass  rate  versus  time,  is  on  the  right.  To  compute  a  mean  residence  time 
based  on  a  moment  analysis  of  the  breakthrough  curve  at  the  extraction  well,  the  following  tem¬ 
poral  moments  were  computed: 

Zero  Absolute  Moment  [M\: 


First  Absolute  Moment  [M  •  t]: 


First  Normalized  Moment  [t]: 


Mo,t 


Computed  Mean  Residence  Time: 


t' 

-^res 


J- pulse 

2 


Where  TpuUe  is  the  duration  of  a  reetangular  tracer  input  pulse.  Tres  and  must  be 
compared  to  determine  whether  our  computational  model  is  self-consistent.  Figure  17  shows 
the  breakthrough  curve  at  the  extraction  wells  resulting  from  a  Tp^ise  =  200  hour  injection 
pulse  of  unit  concentration  and  0.192  ^  injection  rate.  Note  that  the  computation  was  per¬ 
formed  for  7000  hrs  after  which  an  exponential  tail  (dotted  portion  of  the  curve)  was  fit  to  the 
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per  Ubit  Time 


Figure  16.  Schematic  of  Injection  Pulse  and  Mass  Recovery  Curves 
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computed  portion  of  the  breakthrough  curve  as  is  commonly  done  [20,  pg  4]  .  Essentially,  the 
computed  breakthrough  curve  is  plotted  on  semi-Log  coordinates.  The  trailing  region  of  this 
curve  (roughly  the  last  1000  hours)  looks  like  a  straight  line  in  this  coordinate  frame.  An  ex- 
ponetial  tail  (which  also  looks  like  a  straight  line  in  these  coordinates)  can  be  matehed  directly. 
Figure  18  shows  the  same  breakthrough  eurve  with  its  parent  rectangular  injection  pulse.  Log- 
Log  coordinates  are  necessary  for  clarity.  The  partial  overlap  in  the  injection  and  break¬ 
through  profiles  indicates  that  contaminant  had  begun  to  arrive  at  the  extraction  wells  before 
the  completion  of  the  input  pulse.  The  area  under  the  breakthrough  curve  was  calculated  using 
the  trapezoidal  rule  for  the  numerically-generated  portion  of  the  curve  whereas  the  area  under 
the  exponential  tail  was  computed  analytically.  The  computed  result  of  T'^es  =  33834  hours 
compares  very  favorably  with  the  theoretical  Trea  =29161  hours. 

Based  on  the  favorable  result  of  the  moment-analysis  cross-check,  we  can  conclude  that 
our  numerical  formulation  constitutes  a  realistic  model  of  the  injection-extraction  dynamics  in 
the  test  cell  under  investigation. 

In  performing  the  moment  analysis,  cross-sectional  contour  profiles  of  injected  contam¬ 
inant  in  the  test  cell  were  generated  as  a  by-product.  These  are  shown,  for  various  times,  in 
figures  19  through  22.  Colors  represent  aqueous  contaminant  concentration.  Examining  these 
contour  plots  provides  valuable  insight  into  the  nature  of  the  aquifer/aquitard  system.  The  first 
figure  (19)  shows  the  contaminant  plumes  after  the  200  hr  injection  pulse.  Clearly,  the  contami¬ 
nant  has  made  significant  inroads  into  the  aquifer  material  but  little  progress  in  the  aquitard  lay¬ 
ers.  Moving  particularly  rapidly  is  the  contaminant  in  the  orange  medium/coarse  sand  where  the 
flow  rate  is  highest  (figure  15).  At  the  end  of  the  200  hr  contaminant  injection  pulse,  the  clean 
water  flushing  commences  and  in  the  subsequent  plots  we  note  the  lateral  dispersive  spread 
of  contaminant  into  the  aquitard  layers,  even  as  the  flushing  acts  to  force  contaminant  in  the 
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LOG[  Mass  Recovery  Rate  (KigAiir)  ] 


Figure  18.  Injection  Pulse  and  Extrapolated  Extraction  Breakthrough 
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aquifer  layers  towards  the  extraction  wells.  Recall  that  high  hydraulic  conductivity  is  unimpor¬ 
tant  in  diffusion-dominated  processes  such  as  this  lateral  spreading  into  the  aquitard  material. 
Moreover,  since  it  is  not  easy  to  purge  aquitard  layers  of  contaminant  by  clean  water  flushing, 
it  follows  that  once  contaminant  has  ensconced  itself  in  these  aquitard  layers,  it  will  be  difficult 
to  remove  -  hence  the  ’’diminishing  returns”  or  tailing  we  see  in  the  breakthrough  curve  (figure 
17).  Thus,  the  pulse  injection  simulation  has  provided  graphical  and  computational  evidence 
that  the  contaminated  cell  behaves  dispersively  under  clean  water  flushing  and  resists  efforts  to 
purge  it. 
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Figure  19.  Contaminant  Plumes  After  200  hour  Injection  Pulse  (Aqueous  Concentration) 
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4.3  EFFECT  OF  RETARDATION  FACTOR 


Since  the  moment  analysis  simulation  showed  that  the  dispersive  nature  of  the  soil  system 
in  question  played  an  important  role  in  delaying  contaminant  flush-out,  the  question  naturally 
arose  as  to  the  relative  importance  of  the  two  delaying  effects:  dispersion  and  sorption.  Recall 
that  if  the  sorptive  reaction  behaves  linearly  and  is  at  equilibrium,  the  solute  will  move  at  an  av¬ 
erage  velocity  equal  to  the  average  linear  velocity  of  the  ground  water  divided  by  the  retardation 
factor.  Therefore,  retardation  factor  (R)  is  key  to  describing  the  effect  sorption  has  on  the  test 
cell’s  hydrodynamic  characteristics.  The  retardation  factor  for  a  given  soil  type  is  defined  as 


i?=  1  + 


PB  ' 
Ow 


where  pg  is  the  bulk  density  -  the  mass  of  the  solids  divided  by  the  volume  of  the  media.  Kd  is 
the  distribution  coefficient  and  6w  is  the  water  capacity.  The  no-sorption  (i?  =  1)  case  forms 
a  logical  basis  for  comparison  with  the  actual  Dover  simulation.  The  retardation  factor  can  be 
forced  to  one  by  setting  the  distribution  coefficient  Kd  in  turn  equal  to  zero.  Thus,  another 
transport  solution  was  run  for  a  200  hour  injection  pulse,  as  in  section  4.2  above,  but  with 
retardation  factor  equal  to  zero  for  all  soil  layers.  This  was  to  examine  the  effect  sorption  had 
on  contaminant  transport.  Figures  23  and  24  show  the  effect  retardation  factor  has  on  the  mass 
recovery  characteristics  of  the  test  cell.  Namely,  for  i?  =  1  the  breakthrough  curve  peaks  occurs 
28%  sooner  and  peaks  21%  higher.  This  result  was  consistent  with  the  expected  behaviour  for 
this  case,  and  provided  additional  validation  of  the  simulation. 


57 


Mass  Recowjy  Rate  (Kg,4iir) 


4.4  REMEDIATION  BY  CONTINUOUS  PUMPING 


A  sample  input  file  for  the  transport  solution  is  given  in  appendix  B.2.  This  simulation 
modeled  contaminant  purge  when  the  initial  contaminant  distribution  shown  in  figure  25  was 
subjected  to  clean  water  flushing.  This  initial  contaminant  distribution  is  shown  in  profile  in 
figure  26.  The  four  monitoring  points  (refer  to  table  4)  are  also  shown  in  this  plot.  Figure  27 
shows  how  the  contaminant  concentrations  at  the  four  monitoring  points  identified  in  table  4  and 
figure  3  vary  with  time.  There  are  several  features  of  this  plot  that  stand  out.  First  of  all,  there 
is  an  initial  ’’lateral  reconciliation”  of  contaminant  in  which  the  four  curves  converge  rapidly. 
This  indicates  a  smoothing  of  the  peaks  in  the  initial  contaminant  profile  (figure  25)  and  this 
initial  phase  is  over  by  T  =  1 500  hours.  Second  comes  the  ”transition”  phase  of  the  flush.  Here, 
the  ensemble  slope  of  the  reconciled  curves  shallows  out,  implying  a  decreasing  contaminant 
extraction  rate.  In  the  third  ’’mature”  phase,  the  concentration  at  the  former  ’’peak”  (monitor 
point  #4)  falls  below  that  of  the  clay  layer  edge  (monitor  point  #1)  so  that  the  clay  layer  now 
harbors  a  higher  concentration  of  contaminant  than  the  layers  above  it.  This  in  turn  suggests  that 
the  clay  aquitard  stubbornly  retains  contaminant  then  metes  is  out  gradually  -  recontaminating 
the  aquifer  layer  above  and  partially  counteracting  the  convective  purge  occurring  there.  Figure 
28  provides  another  view  of  these  phenomena  and  implies  that  at  high  T  (hours),  there  may  be 
more  of  a  redistribution  of  contaminant  than  a  flush-out. 

Turning  our  attention  now  to  a  time-ordered  series  of  cross-sectional  contour  profiles,  fig¬ 
ures  29  through  36  show  that  the  flushing  action  has  the  effect  of  removing  contaminant  from 
the  cell,  smoothing  down  the  initial  concentration  peaks,  and  laterally  dispersing  contaminant 
into  the  upper  (aquifer)  layers  of  the  cell.  This  is  illustrated  clearly  in  the  line-plot  figure  28. 
Note  also  that  there  is  substantial  qualitative  difference  between  the  situations  depicted  in  the 
cross-sectionals  for  T  =  0  hours,  T  =  200  hours,  and  T  =  400  hours,  but  very  little  qualitative  dif- 
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Figure  25.  Initial  Contaminant  Profile  T  =  0  hrs 
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AQtJEOUS  CONCENTRATION  (Kg /m^S) 

Figure  26.  Inital  Contaminant  Profile  and  4  Monitoring  Points 
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AQUEOUS  CONCENTRATION  (Kg/m®) 


Figure  27.  Contaminant  Concentration  vs  Time  (Hours)  at  the  Four  Monitoring  Points 
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Figure  28.  Time  and  Depth  Variation  of  Contaminant  Concentration  at  Center  of  Cell 
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ference  between  the  profiles  for  T  =  1 200  hours  and  T  =  1400  hours.  This  is  commensurate  with 
the  "tailing”  phenomenon  shown  in  figure  27.  On  another  level,  there  are  several  time  scales 
in  concert  during  the  flushing  dynamic.  To  name  a  few;  a  diffusive  time  scale  inversely  propor¬ 
tional  to  the  fluid  diffusivity  Dm,  a  convective  time  scale  associated  with  the  mean  longitudinal 
velocity  through  the  cell,  and  a  time  scale  based  on  the  mixing  shear  velocity  at  the  interface 
between  the  aquitard  clay  layer  and  it’s  adjacent  aquifer.  These  multiple  time  scales  reflect  the 
complicated  nature  of  the  physics.  They  also  make  the  numerical  simulation  difficult  (stiff)  due 
to  the  range  of  time  resolution  required  for  various  time  steps.  From  an  operational  standpoint, 
the  only  real  conclusion  to  be  drawn  is  that  in  the  presence  of  a  contaminant-retaining  aquitard, 
continuous  pumping  offers  diminishing  returns. 


65 


Figure  31.  Contaminant  Profile 


600  Hours  of  Clean-Water  Flushing 
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Figure  34.  Contaminant  Profile  after  1200  Hours  of  Clean- Water  Flushing 
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Figure  35.  Contaminant  Profile  afterl400  Hours  of  Clean- Water  Flushing 
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Figure  36.  Contaminant  Profile  after  3600  Hours  of  Clean-Water  Flushing 
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4.5  REBOUND 


This  case  examined  what  happens  when  the  injection  pumps  are  switched  off  after  the 
flushing  operation  has  reached  it’s  "mature”  stage.  The  question  was  whether  the  clay  aquitard 
would  continue  to  mete  out  contaminant  diffusively  in  the  absence  of  an  injection-induced  ve¬ 
locity  field.  However,  examining  figure  39,  it  is  evident  that  the  peak  in  the  T=3600  hours 
contaminant  concentration  profile  occurs  between  monitor  points  1  &  2  and  not  within  the  clay 
layer.  Beginning  with  a  normal  pump-down  as  before,  the  pumps  were  switched  off  at  3600 
hours,  allowing  the  plant  to  stand  down  for  19O0  hours.  The  contaminant  level  at  monitoring 
point  #2  (just  above  the  clay  layer)  was  carefully  observed.  As  depicted  in  figure  37,  the  con¬ 
taminant  begins  to  rise  again  before  leveling  out  at  approximately  5500  hours.  This  is  rebound. 
However,  figure  38  shows  that  even  during  pump  stand-down,  the  contaminant  level  at  the  edge 
of  the  clay  layer  (monitoring  point  #1 )  continued  to  fall  at  almost  the  same  rate  as  during  pump- 
on.  This  is  because  MPl  is  affected  by  a  concentration  gradient  favoring  local  decrease.  The 
reason  why  the  concentrations  registering  at  monitor  points  1  &  2  diverge  from  each  other  is  that 
the  concentration  gradients  affecting  these  two  points  are  in  opposite  directions.  Moreover,  the 
difference  in  slope  magnitude  between  these  two  curves  mirrors  the  difference  in  the  two  gradi¬ 
ents  affecting  MPl  and  MP2.  Figure  39  shows  that  MPl  is  at  a  much  higher  concentration  than 
the  point  immediately  beneath  it  in  the  clay  layer  whereas  MP2  is  at  an  only  slightly  lower  con- 
centation  than  the  point  adjacent  to  it.  The  important  point  is  that  mass  transfer  is  by  diffusion. 
Figure  37  shows  that  at  pump-on  at  5500  hours  the  slope  of  the  concentration  vs.  time  curve 
for  monitoring  point  #2  literally  plunges,  with  a  slope  much  steeper  than  that  for  the  portion 
previous  to  pump  stand-down  -  corroborating  the  fact  that  a  quantity  of  mass  arose  diffusively 
from  below  during  stand-down,  availing  itself  for  subsequent  purge  under  the  convective  ac¬ 
tion  of  pump-on.  The  steeper  slope  of  the  second  pump-on  cycle  illustrates  the  increased  mass 
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removal  efficiency  that  the  pulse  pumping  advocates  suggest.  Finally,  note  the  rapid  so-called 
”lateral  reconciliation”  between  contaminant  levels  at  monitoring  points  1  and  2.  This  is  due 
to  the  powerful  shear-induced  convective  mixing  which  results  from  the  steep  velocity  gradient 
between  these  two  relatively  close  layers  (see  figure  1 5  and  table  3). 
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ELEVATION  $n) 


Figure  39.  Contaminant  Concentration  in  the  Vicinity  of  the  4  Monitoring  Points  at  T=3600  Hours 
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4.6  COMPARISON  WITH  FIELD  EXPERIMENTAL  DATA 

Details  of  the  experimental  work  conducted  at  Dover  AFB  are  provided  in  [27]  .  Figure 
Ho  below  is  the  PCE  portion  of  figure  G-1  (Selected  Elution  Curves...)  on  page  188  of  that 
work.  The  accompanying  pump  rate  schedule  does  not  appear  in  the  report  by  Mackay,  Ball  et 
al  [27]  but  was  provided  to  this  author  by  one  of  the  collaborators,  D.M.  Mackay  [26]  .  Figure 
41  below  shows  the  pump  schedule  at  the  field  site.  The  reasons  for  the  variance  in  the  pump 
rate,  according  to  the  experimenters,  included  pump  component  wear-outs,  clogged  filters  and 
valves,  and  a  myriad  of  other  unforseeable  circumstances  that  conspired  against  ideality.  Note 
that  the  average  pump  rate  over  the  entire  3700  hours  shown  is  0.196  ^  which  is  2%  higher 
than  the  0.192  ^  assumed  for  the  present  numerical  work.  Also  note  that  the  experimenters 
recorded  elution  data  points  at  different  times  than  the  pump  strength  data  points  so  it  was 
necessary  to  perform  an  interpolation  between  the  latter  to  provide  pump  rate  data  corresponding 
temporally  to  the  elution  data.  The  inner  product  of  elution  and  pump  strength  data  provided 
the  breakthrough  curve  shown  in  figure  42.  Shown  in  the  same  figure  is  the  breakthrough  curve 
provided  by  FEMWATER.  The  breakthrough  curve  resulting  from  the  numerical  simulation  is 
much  flatter  than  the  experimental  result,  undershooting  at  low  hours  and  overshooting  at  high 
hours.  The  main  reason  for  the  discrepancy  appears  to  be  a  FEMWATER  tendency  to  be  over- 
dispersive.  Clearly,  the  breakthrough  curve  from  the  experiment  shows  plug-like  elution  of 
mass  from  the  test  cell.  FEMWATER’s  solution  on  the  other  hand  seems  to  be  smeared,  with 
dispersion  playing  a  larger  role  than  advection. 
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AQUEOUS  COHCEHTRATIOK(Mraosrmas/l.) 


Figure  40.  Elution  Data  for  PCE  (Fig  G1  in  [27]) 
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Figure  42.  Breakthrough  Curves  for  FEM\\M"ER  and  Field  Experiment 


4.7  COMPARISON  WITH  PREVIOUS  NUMERICAL  DATA 


Herman  [21]  conducted  numerical  simulations  based  on  the  aqueous  concentration  profile 
shown  in  figure  43.  Figure  44  shows  a  comparison  between  Herman’s  and  FEMW^TER  spatial 
concentration  profiles  at  a  point  corresponding  to  ’’monitoring  point  #2”  (table  3).  Again,  Her¬ 
man’s  results  show  a  plug-like  elution  dominated  by  advection  whereas  FEMWATER’s  solution 
contains  smearing  and  mixing  -  hence  the  flatter  breakthrough  curve. 
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ELEVATION  (Metjeis) 


Figure  43.  Initial  Concentration  Profile  from  Herman 
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AQUEOUS  CONCENTRATION  (MicroOr/L) 


Figure  44.  Comparison  of  Spatial  Concentration  Traces 
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4.8  GENERAL  OBSERVATIONS  ON  FEMWATER  PERFORMANCE 


In  the  course  of  performing  the  calculations  described  in  the  previous  sections,  several 
trends  were  noticed: 

1.  The  use  of  the  full  curves  for  soil  hydraulic  properties  (section  3.2.e)  appeared  to  be 
unnecessary:  For  the  present  calculations,  the  test  cell  was  saturated  and  no  benefit  was 
derived  by  representing  the  curves  in  their  unsaturated  ranges.  The  computation  converged 
to  the  identical  solution  with  the  constant  saturated  value  of  each  property  extended  into 
the  unsaturated  portion  of  h. 

2.  Though  the  FEMW\TER  manual  claims  that  a  convergence  tolerance  of  10“^  should 
suffice  for  the  flow  calculation,  the  author  found  that  decreasing  this  to  10~®  affected  the 
magnitude  of  the  computed  velocity  lt;|  as  listed  in  table  5  below.  Clearly  a  24%  increase 
in  mean  |r;|  is  not  inconsequential,  particularly  in  light  of  the  fact  that  the  velocity  field 
plays  such  a  key  role  in  the  transport  equation.  The  rate  of  transport  of  contaminant  is 
sensitive  to  the  magnitude  of  the  velocity.  Among  experienced  researchers  in  the  field, 
Zakikhani  [41]  recounts  at  least  one  instance  in  which  tightening  the  convergence  criterion 
on  velocity  made  the  difference  (in  the  final  transport  solution)  between  poor  agreement 
with  experimental  results  and  excellent  agreement  therewith. 

3.  Tightening  the  convergence  criterion  on  transport  calculations  did  not  seem  to  make  any 
difference  in  computed  transport  solutions  so  the  manual-recommended  value  of  10“^  was 
used  throughout. 


Table  5.  Effects  of  Tightening  Convergence  Tolerance  on  Flow  Solution 


Tole=  10"^ 

Tole=  10-° 

Percent  Increase 

Min  |x;| 

1.4x10"^ 

2.12x10-^ 

51  % 

Max  uj 

0.07308 

0.07471 

2% 

Mean  |v| 

0.00392 

0.00488 

24% 

Std  Dev  lt>| 

0.00672 

0.00727 

8% 

86 


5.  CONCLUSIONS  &  RECOMMENDATIONS  FOR 


FURTHER  STUDY 


5.1  CONCLUSIONS 

The  author  has  succeeded  in  bringing  the  GMS/FEMWATER  package  on  line.  The  Air  Force  now 
has  a  powerful  tool  for  the  analysis  of  ground-water  hydrodynamics  which  is  central  to  the  en¬ 
vironmental  science.  The  following  numbered  statements  summarize  the  conclusions  drawn  from 
the  research  conducted  with  this  tool; 

1.  FEMWATER  is  capable  of  modeling  tailing  and  rebound  without  an  additional  microscale 
model.  Specifically,  it  is  capable  of  representing  diffusion  processes  when  the  soil  reflects 
macroscopic  structure  without  microscale  parameters  which  are  difficult  to  characterize. 

2.  FEMWATER  appears  to  impose  a  lot  of  numerical  dispersion.  As  a  result,  calculations  of 
contaminant  transport  are  dominated  by  smearing  and  mixing  instead  of  advection. 

3.  FEMWATER  has  displayed  an  insensitivity  to  soil  hydraulic  properties  when  the  plant  is 
saturated. 

4.  FEMWATER  has  illustrated  the  effect  of  retardation  factor  on  the  sorptive  retention  of 
contaminant  in  the  test  cell:  Retardation  results  in  a  temporally  delayed  peaking  of  the 
breakthrough  curve  and  a  lower  peak. 

While  the  FEMWATER  results  have  been  shown,  through  the  moment  analysis,  to  be  self- 
consistent,  the  question  of  whether  they  adequately  represent  reality  is  still  open. 

5.2  RECOMMENDATIONS 

Future  work  should  focus  on: 

1.  Understanding  the  reasons  behind  FEMWATER’s  tendency  to  be  over-dispersive. 

2.  More  careful  and  detailed  comparisons  with  actual  field  data.  More  effort  must  be  made  to 
model  the  soil  strata  as  they  really  are,  pursuant  to  better  numerical  simulations. 

3.  Effort  to  improve  the  well  discretization  model  by  doing  away  with  the  assumption  that 
each  source’s  hydraulic  zone  of  influence  extends  to  the  centerline  of  the  adjacent  layer  of 
elements. 
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APPENDIX  A  UNIQUENESS  IN  THE  ONE-DIMENSIONAL 

DIFFUSION  EQUATION 


Consider  the  problem: 


dC 

dt 


d'^C 

dx'^ 


+  F{x,  i) 


C{x,0)  =  f{x) 


=  Oatx  —  0,x  =  L 
ox 


Suppose  now  that  f{x)  =  0  then  our  problem  becomes 


dt 


d^c 

dx^ 


F{x,t) 


C{x,0)  =  0 


=  0  at  X  =  0,  X  =  L 
dx 

For  which  we  recognize  that  if  F{x^  t)  =  0,  then  C  =  constant  is  a  solution  to  the  partial 


differential  equation.  The  implications  of  this  (see  for  example  [8]  )  are  twofold: 

1 .  A  solution  may  exist  for  some  problems  and  not  for  others 

2.  If  a  solution  does  exist  it  may  not  be  unique;  Another  solution  may  exist  within  an  additive 
constant 


For  uniqueness  in  the  above  instance,  there  must  be  a  compatibility  relation  between  the 
forcing  function  and  the  initial/boundary  conditions. 
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APPENDIX  B  SAMPLE  FEMWATER  INPUT  FILES 


The  following  FEMWATER  input  files  are  compatible  with  the  version  1 . 1  of  FEMWATER. 
dated  1  August  1995 


B.1  INPUT  FILE  FOR  FLOW  SOLUTION 

3DFEMWBC 
T1  SMOOTHED  GRID 
T2  TARIQ  HASHIM 
T3  7MAR97 
OPl  10 

OP2  0  0  0  0  1  22 

OP3  l.OOOOOOOOe+00  l.OOOOOOOOe+00  1 .00000006e+00  1  .OOOOOOOOe+00 
OP4  0  0  1 

IPl  40  10  400  l.OOOOOOOOe-06  1  .OOOOOOOOe-06 

IP2  10  10  l.OOOOOOOOe-03  1  .OOOOOOOOe-03 

IP3  10  50  5.00000000e-01  5.00000000e-01 

IP4  l.OOOOOOOOe-03  l.OOOOOOOOe-01 

PTl  1  1  1 2 

TCI  3.00000000e+02 

TC2  0  1.20000000e+00 

TC3  0 

OCl  0  0  0  1 

OC2  0  0 

OC3  1  0  1 

OC4  5  1  2  3  4  5 

MPl  0 

MP3  l.OOOOOOOOe+03  4.68000000e+00  1 .27000000e+08 

MP4  1  .OOOOOOOe+00  O.OOOOOOOe+00  O.OOOOOOOe+00  O.OOOOOOOe+00  1  .OOOOOOOe+00  O.OOOOOOOe+00  O.OOOOOOOe+00 
O.OOOOOOOe+00 

MP2  7  3.60000000e-05  3.60000000e-05  3.60000000e-05  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  7  O.OOOOOOOe+00 1.3000000e+03  2.3000000e-01  2.3000000e-02  2.9800000e-066.0000000e-01  O.OOOOOOOe+00 
l.OOOOOOOe+00 
SPI  7  1  23 

MP2  2  4.32000000e-02  4.32000000e-02  4.32000000e-02  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  2  O.OOOOOOOe+00 1.7000000e+03  2.3000000e-0I  2.3000000e-02  2.9800000e-066.0000000e-01  O.OOOOOOOe+00 
l.OOOOOOOe+00 
SPI 2  1  2  3 

MP2  3  l,08000000e-03  1 .08000000e-03  1 .08000000e-03  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  3  O.OOOOOOOe+00  1.7000000e+03  2.3000000e-01  2.3000000e-022.9800000e-066.0000000e-01  O.OOOOOOOe+00 
l.OOOOOOOe+00 
SPI  3  1  23 

MP2  5  8.28000000e-02  8.28000000e-02  8.28000000e-02  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  5  O.OOOOOOOe+00  1 .7000000e+03  2.3000000e-01  2.3000000e-02  2.9800000e-06  6.0000000e-0 1  O.OOOOOOOe+00 
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l.OOOOOOOe+00 
SPl  5  1  23 

MP2  4  5,40000000e-02  5.40000000e-02  5.40000000e-02  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  4  O.OOOOOOOe+00  1 .7000000e+03  2.3000000e-0 1  2.3000000e-02  2.9800000e-06  6.0000000e-0 1  O.OOOOOOOe+00 
l.OOOOOOOe+00 
SPl 4  I  2  3 

MP2  6  3.60000000e-06  3.60000000e-06  3.60000000e-06  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  6  O.OOOOOOOe+00  1.3000000e+03  2.3000000e-01  2.3000000e-02  2.9800000e-06  6.0000000e-01  O.OOOOOOOe+00 
l.OOOOOOOe+00 
SPl 6  1  2  3 

MP2  8  l.OOOOOOOOe-04  1  .OOOOOOOOe-04  l.OOOOOOOOe-04  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  8  O.OOOOOOOe+00  1 .5000000e+03  23000000e-01  2.3000000e-02  2.9800000e-06  6.0000000e-01  O.OOOOOOOe+00 
LOOOOOOOe+00 
SPl  8  1  23 

MP2  1  1.08000000e-03  1 .08000000e-03  L08000000e-03  0,00000000e+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  1  O.OOOOOOOe+00  1.7000000e+03  2.3000000e-01  23000000e-022.9800000e-066.0000000e-01  O.OOOOOOOe+00 
l.OOOOOOOe+00 
SPl  1  1  2  3 

XYl  1  3  0  0  0  O.OOOOOOOOOOOOOOOe+00  moisture 
-20.000  0.36000 
0.000  0.36000 
20.000  0.36000 

XYl  2  3  0  0  0  O.OOOOOOOOOOOOOOOe+00  conductivity 
-20.000  1.00000 
0.000  1.00000 
20.000  1.00000 

XYl  3  3  0  0  0  O.OOOOOOOOOOOOOOOe+00  watercapacity 
-20.000  0.00000 
0.000  0.00000 
20.000  0.00000 

XYl  4  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1.14867000e-03 

XYl  5  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  l.OOOOOOOOe+00 

XYl  6  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -1.14867000e-03 

XYl  7  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  4.25430000e-05 

XYl  8  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -4.25430000e-05 

XYl  9  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1.1 6926  lOOe-03 

XYl  10  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -1.1 6926  lOOe-03 

XYl  11  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  2.26897700e-03 

XYl  12  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
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0.000  -2.26897700e-03 

XYl  13  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  7.74100000e-05 

XYl  14  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -7.741  OOOOOe-05 

XYl  15  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  9.49090000e-05 

XYl  16  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -9.49090000e-05 

XYl  17  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1.16370000e-04 

XYl  18  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -1.16370000e-04 

XYl  19  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1.42681000e-04 

XYl  20  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000-1. 4268  lOOOe-04 

XYl  21  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1.74929000e-04 

XYl  22  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -1.74929000e-04 

XYl  23  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  3.23226220e-03 

XYl  24  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -3.23226200e-03 

XYl  25  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  6.16594500e-03 

XYl  26  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -6.16594500e-03 

XYl  27  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  5.9570  lOOOe-03 

XYl  28  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -5.95701000e-03 

XYl  29  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  5.75516500e-03 

XYl  30  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -5.755 16500e-03 

XYl  31  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  5.56088400e-03 

XYl  32  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -5.56088400e-03 

XYl  33  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  6.78077400e-03 

XYl  34  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -6.78077400e-03 

XYl  35  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  7.09446000e-03 


91 


XYl  36  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -7.09446000e-03 

XYl  37  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  5.3375 1600e-03 

XYl  38  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -5.3375 1600e-03 

XYl  39  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  4.0 1473400e-03 

XYl  40  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -4.01473400e-03 

XYl  41  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  3.02101 700e-03 

XYl  42  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -3.02101700e-03 

XYl  43  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  2.27301  lOOe-03 

XYl  44  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -2.27301100e-03 

XYl  45  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1.70983200e-03 

XYl  46  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -1.70983200e-03 

XYl  47  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1. 28654  lOOe-03 

XYl  48  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000-1. 28654  lOOe-03 

XYl  49  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  5.54708000e-04 

XYl  50  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -5.54708000e-04 

XYl  51  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  2.03920000e-05 

XYl  52  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 

0.000  -2.03920000e-05 

PSl  32147  4 

PS2  32147  5 

PSl  32146  6 

PSl  32145  4 

PS2  32145  5 

PSl  32144  6 

PSl  32143  4 

PS2  32143  5 

PSl  32142  6 

PSl  311107 

PS2311105 

PSl  31109  8 

PSl  31108  7 


PS2  31108  5 
PSl  31107  8 
PSl  31106  7 
PS2  31106  5 
PSl  31105  8 
PSl  30073  9 
PS2  30073  5 
PSl  30072  10 
PSl  30071  9 
PS2  30071  5 
PSl  30070  10 
PSl  30069  9 
PS2  30069  5 
PSl  30068  10 
PSl  29036  11 
PS2  29036  5 
PSl  29035  12 
PSl  29034  11 
PS2  29034  5 
PSl  29033  12 
PSl  29032  11 
PS2  29032  5 
PSl  29031  12 
PSl  27999  13 
PS2  27999  5 
PSl  27998  14 
PSl  27997  13 
PS2  27997  5 
PSl  27996  14 
PSl  27995  13 
PS2  27995  5 
PSl  27994  14 
PSl  26962  15 
PS2  26962  5 
PSl  26961  16 
PSl  26960  15 
PS2  26960  5 
PSl  26959  16 
PSl  26958  15 
PS2  26958  5 
PSl  26957  16 
PSl  25925  17 
PS2  25925  5 
PSl  25924  18 
PSl  25923  17 
PS2  25923  5 
PSl  25922  18 


PSl  25921  17 
PS2  25921  5 
PSl  25920  18 
PSl  24888  19 
PS2  24888  5 
PSl  24887  20 
PSl  24886  19 
PS2  24886  5 
PSl  24885  20 
PSl  24884  19 
PS2  24884  5 
PSl  24883  20 
PSl  23851  21 
PS2  23851  5 
PSl  23850  22 
PSl  23849  21 
PS2  23849  5 
PSl  23848  22 
PSl  23847  21 
PS2  23847  5 
PSl  23846  22 
PSl  22814  23 
PS2  22814  5 
PSl  22813  24 
PSl  22812  23 
PS2  22812  5 
PSl  22811  24 
PSl  22810  23 
PS2  22810  5 
PSl  22809  24 
PSl  21777  25 
PS2  21777  5 
PSl  21776  26 
PSl  21775  25 
PS2  21775  5 
PSl  21774  26 
PSl  21773  25 
PS2  21773  5 
PSl  21772  26 
PSl  20740  27 
PS2  20740  5 
PSl  20739  28 
PSl  20738  27 
PS2  20738  5 
PSl  20737  28 
PSl  20736  27 
PS2  20736  5 


PSl  20735  28 
PSl  19703  29 
PS2  19703  5 
PSl  19702  30 
PSl  19701  29 
PS2  19701  5 
PSl  19700  30 
PSl  19699  29 
PS2  19699  5 
PSl  19698  30 
PSl  18666  31 
PS2  18666  5 
PSl  18665  32 
PSl  18664  31 
PS2  18664  5 
PSl  18663  32 
PSl  18662  31 
PS2  18662  5 
PSl  18661  32 
PSl  17629  33 
PS2  17629  5 
PSl  17628  34 
PSl  17627  33 
PS2  17627  5 
PSl  17626  34 
PSl  17625  33 
PS2  17625  5 
PSl  17624  34 
PSl  16592  35 
PS2  16592  5 
PSl  16591  36 
PSl  16590  35 
PS2  16590  5 
PSl  16589  36 
PSl  16588  35 
PS2  16588  5 
PSl  16587  36 
PSl  15555  37 
PS2  15555  5 
PSl  15554  38 
PSl  15553  37 
PS2  15553  5 
PSl  15552  38 
PSl  15551  37 
PS2  15551  5 
PSl  15550  38 
PSl  14518  39 
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PS2  14518  5 
PSl  14517  40 
PSl  14516  39 
PS2  14516  5 
PSl  14515  40 
PSl  14514  39 
PS2  14514  5 
PSl  14513  40 
PSl  13481  41 
PS2  13481  5 
PSl  13480  42 
PSl  13479  41 
PS2  13479  5 
PSl  13478  42 
PSl  13477  41 
PS2  13477  5 
PSl  13476  42 
PSl  12444  43 
PS2  12444  5 
PSl  12443  44 
PSl  12442  43 
PS2  12442  5 
PSl  12441  44 
PSl  12440  43 
PS2  12440  5 
PSl  12439  44 
PSl  11407  45 
PS2  11407  5 
PSl  11406  46 
PSl  11405  45 
PS2  11405  5 
PSl  11404  46 
PSl  11403  45 
PS2  11403  5 
PSl  11402  46 
PSl  10370  47 
PS2  10370  5 
PSl  10369  48 
PSl  10368  47 
PS2  10368  5 
PSl  10367  48 
PSl  10366  47 
PS2  10366  5 
PSl  10365  48 
PSl  2074  49 
PS2  2074  5 
PSl  2073  50 
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PSl  2072  49 

PS2  2072  5 

PSl  2071  50 

PSl  2070  49 

PS2  2070  5 

PSl  2069  50 

PSl  1037  51 

PS2  1037  5 

PSl  1036  52 

PSl  1035  51 

PS2  1035  5 

PSl  1034  52 

PSl  1033  51 

PS2  1033  5 

PSl  1032  52 

ICHO  l.OOOOOOOOe+01 

ICCO  l.OOOOOOOOe+01 

ICSO 

ICT  O.OOOOOOOOe+00 

ICF  0  0  0 

END 
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B.2  INPUT  FILE  FOR  TRANSPORT  SOLUTION 


3DFEMWBC 
T1  SMOOTHED  GRID 
T2  TARIQ  HASHIM 
T3  7MAR97 
OPl  1 

OP2  1  1  0  0  1  22 

OP3  l.OOOOOOOOe+00  1  .OOOOOOOOe+00  1  .OOOOOOOOe+00  TOOOOOOOOe+00 
OP4  0  0  1 

IPl  40  10  400  l.OOOOOOOOe-06  1  .OOOOOOOOe-06 
IP2  10  10  l.OOOOOOOOe-03  l.OOOOOOOOe-03 
IP3  10  50  5.00000000e-01  5.00000000e-01 
IP4  TOOOOOOOOe-03  l.OOOOOOOOe-01 
PTl  1  1  1 2 
TCI  1400.0 

TC2  0  l.OOOOOOOOe+00 

TC3  0 

OCl  000  1 

OC2  0  0 

OC3  1  0  200.0 

OC4  5  1  2  3  4  5 

MPl  0 

MP3  l.OOOOOOOOe+03  4.68000000e+00  1.27000000e+08 

MP4  1  .OOOOOOOe+00  O.OOOOOOOe+00  O.OOOOOOOe+00  O.OOOOOOOe+00  1  .OOOOOOOe+00  O.OOOOOOOe+00  O.OOOOOOOe+00 
O.OOOOOOOe+00 

MP2  7  3.60000000e.05  3.600000006-05  3.60000000e-05  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  7  O.OOOOOOOe+00  1 .3000000e+03  2.3000000e-01  2.3000000e-02  2.9800000e-06  6.0000000e-01  O.OOOOOOOe+00 
l.OOOOOOOe+00 
SPI  7  1  23 

MP2  2  4.32000000e-02  4.32000000e-02  4.32000000e-02  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  20.0000000e+00  1.7000000e+03  2.3000000e-01  2.3000000e-02  2.9800000e-066.0000000e-01  O.OOOOOOOe+00 
l.OOOOOOOe+00 
SPI  2  1  23 

MP2  3  1.08000000e-03  1.08000000e-03  1.08000000e-03  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  3  O.OOOOOOOe+00  1 .7000000e+03  2.3000000e-01  2.3000000e-02  2.9800000e-06  6.0000000e-01  O.OOOOOOOe+00 
l.OOOOOOOe+00 
SPI  3  1  23 

MP2  5  8.28000000e-02  8.28000000e-02  8.28000000e-02  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  5  O.OOOOOOOe+00  1 .7000000e+03  2.3000000e-01  2.3000000e-02  2.9800000e-06  6.0000000e-01  O.OOOOOOOe+00 
l.OOOOOOOe+00 
SPI  5  1  23 

MP2  4  5,40000000e-02  5.40000000e-02  5.40000000e-02  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  4  O.OOOOOOOe+00  1.7000000e+03  2.3000000e-01  2.3000000e-02  2.9800000e-06  6.0000000e-01  O.OOOOOOOe+00 
l.OOOOOOOe+00 
SPI  4  1  23 

MP2  6  3.60000000e-06  3.60000000e-06  3.60000000e-06  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
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MP5  60.0000000e+00  1.3000000e+03  2.3000000e-01  2.3000000e-022.9800000e-066.0000000e-01  O.OOOOOOOe+00 

l.OOOOOOOe+00 

SPl  6  1  23 

MP2  8  l.OOOOOOOOe-04  1  .OOOOOOOOe-04  LOOOOOOOOe-04  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  80.0000000e+00  1 .5000000e+03  2.3000000e-01  2.3000000e-022.9800000e-066.0000000e-0l  0,0000000e+00 
l.OOOOOOOe+00 
SPl  8  1  23 

MP2  1  1.08000000e-03  1 .08000000e-03  1 .08000000e-03  O.OOOOOOOOe+00  O.OOOOOOOOe+00  O.OOOOOOOOe+00 
MP5  1  O.OOOOOOOe+00  1.7000000e+03  2.3000000e-01  2.3000000e-02  2.9800000e-066.0000000e-01  O.OOOOOOOe+00 
l.OOOOOOOe+00 
SPl  1  1  2  3 

XYl  1  3  0  0  0  O.OOOOOOOOOOOOOOOe+00  moisture 
-20.000  0.36000 
0.000  0.36000 
20.000  0.36000 

XYl  2  3  0  0  0  O.OOOOOOOOOOOOOOOe+00  conductivity 
-20.000  1.00000 
0.000  1.00000 
20.000  1.00000 

XYl  3  3  0  0  0  O.OOOOOOOOOOOOOOOe+00  watercapacity 
-20.000  0.00000 
0.000  0.00000 
20.000  0.00000 

XYl  4  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  I.14867000e-03 

XYl  5  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  l.OOOOOOOOe+00 

XYl  6  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -1.14867000e-03 

XYl  7  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  4.25430000e-05 

XYl  8  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -4.25430000e-05 

XYl  9  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1.1 6926  lOOe-03 

XYl  10  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -1.16926100e-03 

XYl  11  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  2.26897700e-03 

XYl  12  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -2.26897700e-03 

XYl  13  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  7.74100000e-05 

XYl  14  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -7.74100000e-05 

XYl  15  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  9.49090000e-05 
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XYl  16  1  0  0  0  0,000000000000000e+00  constant 
0.000  -9.49090000e-05 

XYl  17  10  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1.16370000e-04 

XYl  18  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -1.1 63  70000e-04 

XYl  19  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1.42681000e-04 

XYl  20  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -1.42681000e-04 

XYl  21  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1.74929000e-04 

XYl  22  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -1.74929000e-04 

XYl  23  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  3.23226220e-03 

XYl  24  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -3.23226200e-03 

XYl  25  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  6.16594500e-03 

XYl  26  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -6.16594500e-03 

XYl  27  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  5.95701000e-03 

XYl  28  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -5.95701000e-03 

XYl  29  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  5.755 16500e-03 

XYl  30  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -5.755 16500e-03 

XYl  31  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  5.56088400e-03 

XYl  32  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -5.56088400e-03 

XYl  33  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  6.78077400e-03 

XYl  34  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -6.78077400e-03 

XYl  35  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  7.09446000e-03 

XYl  36  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -7.09446000e-03 

XYl  37  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  5.33751600e-03 

XYl  38  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -5.3375 1600e-03 

XYl  39  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 


0.000  4.01473400e-03 

XYl  40  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -4.01473400e-03 

XYl  41  I  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  3.02101700e-03 

XYl  42  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -3.02101700e-03 

XYl  43  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  2.27301100e-03 

XYl  44  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -2.27301  lOOe-03 

XYl  45  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1.70983200e-03 

XYl  46  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -1.70983200e-03 

XYl  47  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  1.28654100e-03 

XYl  48  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000 -1.28654100e-03 

XYl  49  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  5.54708000e-04 

XYl  50  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  -5.54708000e-04 

XYl  51  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 
0.000  2.03920000e-05 

XYl  52  1  0  0  0  O.OOOOOOOOOOOOOOOe+00  constant 

0.000  -2.03920000e-05 

PSl  32147  4 

PS2  32147  5 

PSl  32146  6 

PSl  32145  4 

PS2  32145  5 

PSl  32144  6 

PSl  32143  4 

PS2  32143  5 

PSl  32142  6 

PSl  311107 

PS2311105 

PSl  31109  8 

PSl  31108  7 

PS2  31108  5 

PSl  31107  8 

PSl  311067 

PS2  31106  5 

PSl  31105 8 

PSl  30073  9 

PS2  30073  5 


PSl  30072  10 
PSl  30071  9 
PS2 30071  5 
PSl  30070  10 
PSl  30069  9 
PS2  30069  5 
PSl  30068  10 
PSl  29036  11 
PS2  29036  5 
PSl  29035  12 
PSl  29034  11 
PS2  29034  5 
PSl  29033  12 
PSl  29032  11 
PS2  29032  5 
PSl  29031  12 
PSl  27999  13 
PS2  27999  5 
PSl  27998  14 
PSl  27997  13 
PS2  27997  5 
PSl  27996  14 
PSl  27995  13 
PS2  27995  5 
PSl  27994  14 
PSl  26962  15 
PS2  26962  5 
PSl  26961  16 
PSl  26960  15 
PS2  26960  5 
PSl  26959  16 
PSl  26958  15 
PS2  26958  5 
PSl  26957  16 
PSl  25925  17 
PS2  25925  5 
PSl  25924  18 
PSl  25923  17 
PS2  25923  5 
PSl  25922  18 
PSl  25921  17 
PS2  25921  5 
PSl  25920  18 
PSl  24888  19 
PS2  24888  5 
PSl  24887  20 
PSl  24886  19 


PS2  24886  5 
PSl  24885  20 
PSl  24884  19 
PS2  24884  5 
PSl  24883  20 
PSl  23851  21 
PS2  23851  5 
PSl  23850  22 
PSl  23849  21 
PS2  23849  5 
PSl  23848  22 
PSl  23847  21 
PS2  23847  5 
PSl  23846  22 
PSl  22814  23 
PS2  22814  5 
PSl  22813  24 
PSl  22812  23 
PS2  22812  5 
PSl  22811  24 
PSl  22810  23 
PS2  22810  5 
PSl  22809  24 
PSl  21777  25 
PS2  21777  5 
PSl  21776  26 
PSl  21775  25 
PS2  21775  5 
PSl  21774  26 
PSl  21773  25 
PS2  21773  5 
PSl  21772  26 
PSl  20740  27 
PS2  20740  5 
PSl  20739  28 
PSl  20738  27 
PS2  20738  5 
PSl  20737  28 
PSl  20736  27 
PS2  20736  5 
PSl  20735  28 
PSl  19703  29 
PS2  19703  5 
PSl  19702  30 
PSl  19701  29 
PS2  19701  5 
PSl  19700  30 


PSl  19699  29 
PS2  19699  5 
PSl  19698  30 
PSl  18666  31 
PS2  18666  5 
PSl  18665  32 
PSl  18664  31 
PS2  18664  5 
PSl  18663  32 
PSl  18662  31 
PS2  18662  5 
PSl  18661  32 
PSl  17629  33 
PS2  17629  5 
PSl  17628  34 
PSl  17627  33 
PS2  17627  5 
PSl  17626  34 
PSl  17625  33 
PS2  17625  5 
PSl  17624  34 
PSl  16592  35 
PS2  16592  5 
PSl  16591  36 
PSl  16590  35 
PS2  16590  5 
PSl  16589  36 
PSl  16588  35 
PS2  16588  5 
PSl  16587  36 
PSl  15555  37 
PS2  15555  5 
PSl  15554  38 
PSl  15553  37 
PS2  15553  5 
PSl  15552  38 
PSl  15551  37 
PS2  15551  5 
PSl  15550  38 
PSl  14518  39 
PS2  14518  5 
PSl  1451740 
PSl  14516  39 
PS2  14516  5 
PSl  14515  40 
PSl  14514  39 
PS2  14514  5 


PSl  14513  40 
PSI  13481  41 
PS2  13481  5 
PSl  13480  42 
PSl  13479  41 
PS2  13479  5 
PSI  13478  42 
PSl  13477  41 
PS2  13477  5 
PSl  13476  42 
PSl  12444  43 
PS2  12444  5 
PSl  12443  44 
PSl  12442  43 
PS2  12442  5 
PSl  12441  44 
PSl  12440  43 
PS2  12440  5 
PSl  12439  44 
PSl  11407  45 
PS2  11407  5 
PSl  11406  46 
PSl  11405  45 
PS2  11405  5 
PSl  11404  46 
PSl  11403  45 
PS2  11403  5 
PSl  11402  46 
PSl  10370  47 
PS2  10370  5 
PSl  10369  48 
PSl  10368  47 
PS2  10368  5 
PSl  10367  48 
PSI  10366  47 
PS2  10366  5 
PSl  10365  48 
PSl  2074  49 
PS2  2074  5 
PSl  2073  50 
PSl  2072  49 
PS2  2072  5 
PSl  2071  50 
PSl  2070  49 
PS2  2070  5 
PSI  2069  50 
PSl  1037  51 
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PS2  1037  5 

PSl  1036  52 

PSl  1035  51 

PS2  1035  5 

PSl  1034  52 

PSl  1033  51 

PS2  1033  5 

PSl  1032  52 

ICHO  l.OOOOOOOOe+01 

ICC  0  O.OOOOOOOOe+00 

ICSO 

ICT  O.OOOOOOOOe+00 

ICF  0  1  0 

END 
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